Equivalent circuit model of radiative heat transfer
Stanislav I. Maslovski, Constantin R. Simovski, Sergei A. Tretyakov
aa r X i v : . [ phy s i c s . op ti c s ] N ov Equivalent circuit model of radiative heat transfer
Stanislav I. Maslovski, Constantin R. Simovski, and Sergei A. Tretyakov Departamento de Engenharia Electrot´ecnicaInstituto de Telecomunica¸c˜oes, Universidade de CoimbraP´olo II, 3030-290 Coimbra, Portugal Aalto University, School of Electrical EngineeringSMARAD Center of Excellence, P.O. Box 13000, 00076 Aalto, Finland (Dated: August 15, 2018)Here, we develop a theory of radiative heat transfer based on an equivalent electrical networkrepresentation for the hot material slabs in an arbitrary multilayered environment with arbitrarydistribution of temperatures and electromagnetic properties among the layers. Our approach is fullyequivalent to the known theories operating with the fluctuating current density, while being signifi-cantly simpler in analysis and applications. A practical example of the near-infrared heat transferthrough the micron gap filled with an indefinite metamaterial is considered using the suggestedmethod. The giant enhancement of the transferred heat compared to the case of the empty gap isshown.
PACS numbers: 44.40.+a, 78.67.-n, 42.25.Bs
I. INTRODUCTION
There are two most important heat exchange processes known: thermal conductivity , associated with collectiveoscillations of atoms (i.e., phonons) or electrons in solids including metals, or with convection in fluids or gases, and radiative heat transfer which is associated to the electromagnetic radiation produced by thermally agitated atoms,e.g., the black body radiation. In this work we concentrate on the latter process which can be dominant when thebodies that exchange heat are separated by gaps (empty or filled with a heterogeneous material that weakly conductsthe heat).As is known, thermal radiation results from fluctuations of charge and current density in matter. This phenomenonis governed by the fluctuation-dissipation theorem (FDT) that relates the mean-square fluctuations of a physicalquantity to the dissipation associated with the dynamics of the same quantity. Because in dielectrics the loss isrepresented by the imaginary part of the permittivity, which may be in turn related to the effective conductivity ofthe material, the FDT requires the volumetric current density within a dielectric or conducting body to fluctuate.The FDT constitutes the basis of the present day radiative heat transfer theories which deal with the fluctuatingcurrents and treat them as the principal source of thermal radiation. Such a picture places the radiative heat transfercalculations in the framework of classical electromagnetic theory based on the macroscopic Maxwell equations.Indeed, the well-known theory of radiative heat transfer through narrow vacuum gaps by Polder and van Hove belongs to this class. Similar techniques have been recently developed for the cases when radiative heat transferthrough micron and even submicron gaps is assisted with nanostructured metamaterials (see, e.g., in Refs. 3–5). Inthese approaches the fluctuating current density in the two neighboring media is first obtained from the FDT, followingthe methodology introduced by Rytov. Next, the electromagnetic field produced by the fluctuating currents iscalculated, after which the mean value of the power flow (Poynting vector) across the gap (filled with the metamaterial)is found either with the multiple reflection method, or with a more general transfer matrix approach. Historically, however, thermal agitation of fluctuating currents was first discovered by Johnson in electric circuits and the theory of such thermal fluctuations was developed by Nyquist using ideas which were natural for an engineerdealing with networks of lumped elements and transmission lines. Although this theory was a precursor to the FDT,it is still in wide use in the theory of thermal noise at radio and microwave frequencies. The famous Nyquist’s resultstates that, in any linear passive two-pole (i.e., a single port device with an input represented by two electric contacts)operating at a temperature T the electric thermal fluctuations (in other words, the thermal noise) concentrated withina narrow frequency interval ∆ ν can be equivalently represented by the fluctuating electromotive force (EMF) e ( t ),with the mean-square of fluctuations e = 4Θ( ν, T ) R ( ν )∆ ν, (1)where Θ( ν, T ) = hν/ [exp( hν/k B T ) −
1] is Planck’s mean energy of a harmonic oscillator with h and k B being thePlanck and Boltzmann constants, respectively, and R ( ν ) is the input resistance (real part of the input impedance) ofthe two-pole. The equivalent EMF is then understood as connected in series with the two-pole, which can be nowconsidered noiseless.The beauty of this result is in that no knowledge of the internal structure of the electric network is required, andthat all the information is contained within just a single parameter: the real part of the frequency-dependent inputimpedance. In the terminology of FDT, the equivalent EMF in Nyquist’s formula has the role of a generalized force,and the two-pole input impedance is related to the generalized susceptibility of the system. In other words, theNyquist formula can be obtained by a direct application of FDT to the electric circuit. Note that the fluctuatingcurrent appears in this model only as a reaction to a finite number of lumped voltage sources.In contrast, in the thermal transfer theory based on the full-wave electromagnetic formulation through the fluctu-ating current density understood as a volume-distributed source, the internal structure of the interacting bodies hasto be considered during most of the calculations. Thus, a thermal transfer problem appears in this formulation as aproblem with infinite number of degrees of freedom. The final result, however, happens to be expressed in quantitiesthat abstract away the internal structure, such as the reflection coefficients of material half-spaces in the theory ofPolder and van Hove. This observation suggests that introducing the distributed fluctuating current can be avoidedin many cases. For example, in this work we prove that in stratified media the input impedance concept and theoriginal Nyquist theory can be generalized and used not only in problems related to electric networks, but also infull-wave radiative heat transfer problems. This allows for a significant reduction in complexity of the analysis andmakes the theory readily available for practical calculations.It has to be mentioned that an expression for the equivalent lumped EMF of the thermal noise for the case of alossy material half-space was first derived by Rytov using an approach based on the fluctuating current density. Thisresult constitutes the basis of the thermal noise theory of aperture antennas. Rytov also worked on equivalentfour-pole network representation of a hot material slab. Why then these results have not been widely used in theheat transfer problems? Perhaps, it is because the concepts of input impedance and the equivalent circuit descriptionfor full-wave electromagnetic problems are largely unknown among theorists working in the field of radiative heattransfer. In applied electromagnetics, however, it is well-known that stratified media can be very efficiently treatedwithin the so-called vector transmission line theory (VTLT) which, in essence, assigns an equivalent transmissionline network to every electromagnetic mode (propagating or evanescent) in the system. We would like to stress herethat the VTLT is not an approximation: it is a direct consequence of the Maxwell equations when modal expansionis applied to the electromagnetic field in layered structures. The VTLT allows also for a systematic treatment ofuniaxial and bi-anisotropic media.In this work we extend the VTLT in order to include the effect of the fluctuating current density within the layers.In contrast to a few numerical discrete-element models currently available from the literature (e.g., Ref. 12), thetheory that we develop here is fully analytical. The generalized theory allows us to prove a complete equivalencebetween a volumetric multilayered structure and its circuit theory counterpart, which may be visualized as a chain oftransmission line segments with equivalent fluctuating voltage sources connected at the ports. When concerned withthe radiative heat transfer between the layers, we show that this equivalent network may be reduced to just a seriesconnection of a number of voltage sources representing the fluctuating EMFs and equivalent impedances (each can beunder different temperature), thus, recovering in this way the famous Nyquist result, generalized here to the full-waveelectromagnetic processes in stratified media. Therefore, the calculation of the radiative heat transfer between thelayers reduces in our theory to a number of equivalent circuit theory calculations, which are relatively simple and verysimilar to what is typically done when considering thermal noise in practical electric networks. II. IMPEDANCE REPRESENTATION OF THE CLASSICAL HEAT TRANSFER FORMULA
In this section we establish a connection between the classic Polder–van Hove theory of radiative thermal transferand its representation in terms of the input impedances of the material half-spaces. Formulas by Polder and van Hovefor the density of the radiative heat flux (i.e., power flux of the thermal radiation) across the gap between two thickdielectric slabs [the geometry is defined in Fig. 1(a); the slabs are approximated by half-spaces] read as (our notationscorrespond to Refs. 16–18): S t = ∞ Z dω [Θ( ω, T ) − Θ( ω, T )] M, Θ( ω, T i ) = ~ ωe ~ ωk B Ti − . (2)Here ~ = h/ (2 π ), and T i is the absolute temperature of the i -th medium ( i = 1 , i = 2 represents the gap). If T > T the total radiative flux is directed from medium 1 to medium 3 and can we written as S t = S → − S → where S → is the heat flux produced by medium 1 and absorbed in medium 3 and S → is the flux produced by medium 3and absorbed in medium 1. M from (2) is called the radiative heat transfer function. It depends on the opticalproperties of all three media and may be written as the sum M = M p + M e , where M p and M e are contributions of (a) (b) FIG. 1: (Color online) (a) – Illustration to the general problem formulation. (b) – A possible implementation of medium2 suggested in Ref. 5 transforms the gap into a layered structure filled with so-called hyperbolic metamaterial (see, e.g., inRefs. 19–22) formed by carbon nanotubes or metal nanowires. Interdigital arrangement of nanotubes (nanowires) helps keepradiative heat transfer dominating over the thermal conductance through the gap. propagating and evanescent waves, respectively [note that in this paper we assume the time dependence exp(+ jωt )with j = √− M p = 1 π k Z N p ( ω, q ) q dq, N p ( ω, q ) = (1 − | Γ ( q, ω ) | )(1 − | Γ ( q, ω ) | )4 | − e − jβd Γ ( q, ω )Γ ( q, ω ) | , (3) M e = 1 π ∞ Z k N e ( ω, q ) q dq, N e ( ω, q ) = Im[Γ ( q, ω )]Im[Γ ( q, ω )] e − | β | d | − e − | β | d Γ ( q, ω )Γ ( q, ω ) | , (4)where Γ and Γ are the reflection coefficients of a plane wave harmonic having the spatial frequency (transverse wavenumber) q ≡ | k t | and being incident from a lossless medium 2 (here, medium 2 is free space) to the surfaces of media1 and 3, respectively, β = p k − q is the normal component of the wave vector in medium 2, and k = ω √ ε µ isthe wave number in medium 2. Function N ( ω, q ) is called the spatial spectrum of the radiative heat transfer function.This function was studied in Ref. 16 for the case when photon tunnelling through the vacuum gap was enhancedby surface plasmon-polaritons excited at the gap boundaries. It was shown that the absolute maximum of N ( ω, q )achievable at certain values of ω and q equals 1 / and Γ in (3) and (4) through wave impedances Z , , of the media inregions 1, 2 and 3. The wave impedance defines the ratio between the transverse electric and magnetic fields in aplane wave of a given polarization, i.e., in a given harmonic of the spatial spectrum (see, e.g., in Ref. 11). In termsof the wave impedances, Γ i = Z i − Z Z i + Z , i = 1 , . (5)For isotropic dielectrics the wave impedances are given by the following expressions (see, e.g., in Ref. 11) for TM-wavesand TE-waves, respectively: Z TM i = η β i k ε i , Z TE i = η k β i , (6)where η is the characteristic impedance of free space η = p µ /ε , and β i denotes the normal component of thewave vector in the i -th medium: β i = p k ε i − q . If medium 2 is free space (as it is assumed in the classical theoryof Ref. 2) ε i = 1. For propagating waves q < k and Z is real, for evanescent waves q > k and Z is imaginary. It isimportant that the impedance representation (5) of the reflection coefficients is general for all spatial frequencies q .Substituting Eq. (5) into Eq. (3) we obtain after a rather simple algebra: N p = 16 R R R | ( Z + Z )( Z + Z ) + ( Z − Z )( Z − Z ) e − jβ d | . (7)Here and below we denote R i ≡ Re( Z i ) and X i ≡ Im( Z i ). Substituting Eq. (5) into Eq. (4) we easily deduce: N e = 4 X R R e − | β | d | ( Z + Z )( Z + Z ) + ( Z − Z )( Z − Z ) e − | β | d | . (8)In this case β has imaginary value and it is taken into account that Z = jX . Since X = 0 for propagating wavesand R = 0 for evanescent waves, Eqs. (7) and (8) can be unified into an expression suitable for both regions q < k and q > k : N ( ω, q ) = 4 | Z | R R | e − jβ d | | ( Z + Z )( Z + Z ) + ( Z − Z )( Z − Z ) e − jβ d | . (9)The heat flux density originated from the medium 1 and absorbed in the medium 3 can now be represented asintegral over the complete spatial spectrum (both propagating and evanescent): S → = 1 π ∞ Z dω ∞ Z qdq Θ( ω, T ) N ( ω, q ) . (10)As we show in the next section, Eq. (9) for the spatial spectrum of the heat transfer function N ( ω, q ) can be derivedfrom the VTLT in a way that allows for a straightforward generalization of the results of Polder and van Hove to thecase of stratified uniaxial magneto-dielectric media without a need to introduce the distributed fluctuating currents. III. RADIATIVE HEAT TRANSFER RESULTING FROM AN EQUIVALENT CIRCUIT APPROACH
The possibility to express the radiative heat transfer through the wave impedances of the material layers providesus with an evidence that an equivalent circuit model can be formulated for this problem. Such a model is derivedrigorously in Appendix directly from the Maxwell equations, which results in the VTLT generalized to include theeffect of thermal fluctuations. Throughout this section, however, we will use simple physical reasoning when possible,in order to keep the derivations easy to grasp.As in the previous section, we decompose the fluctuating electromagnetic field into plane waves characterized witha certain polarization state, angular frequency ω , and the transverse wave vector k t . For generality, let us assumethat all materials taking part in the heat transfer are optically uniaxial magneto-dielectric media. It is known thatin uniaxial magneto-dielectrics the independent polarization states correspond to the TE and TM plane waves. Thesame holds for a multilayered structure composed of uniaxial magneto-dielectric layers (some of them can be isotropicor even vacuum gaps) under the condition that the anisotropy axes of all layers are aligned. In such a structure thetwo polarizations are completely independent and can be considered separately. In the following we assume that theanisotropy axis of the layers coincides with the z -axis which is perpendicular to the layers.In Appendix we prove that, in a given layer of the considered multilayered structure ( i -th layer) being under thetemperature T i the transverse components of the time-harmonic fluctuating electric and magnetic fields at the layerinterfaces (labeled here with subscripts 1 and 2) are related as follows Z i Z i Z i Z i ! · (cid:18) n × H i t n × H i t (cid:19) − (cid:18) E i t E i t (cid:19) = 1 √ A (cid:18) e i e i (cid:19) (11)where Z imn = Z imn ( ω, k t ) are the dyadic Z -parameters of the chosen layer, and e i , are the vectorial fluctuating EMFsequivalently representing the thermal-electromagnetic fluctuations within the same layer, and n , are the externalunit normals at the interfaces of the layer. The meaninig of the factor 1 / √ A is explained further in the text. Eq. (11)generalizes the known result of the VTLT to the case of non-vanishing thermal fluctuations. When the right-handside of (11) vanishes this equation represents the definition of the impedance matrix of a passive material layer.The result (11) is obtained in a dyadic form and, thus, is applicable to any polarization state of the electromagneticfield in an anisotropic (not only uniaxial) layer. However, when the states split into the TE and TM waves it ismore convenient to work with scalar Z -parameters — elements of the 2 × ε = ε ⊥ i I t + ε k i z z and the permeability dyadic µ = µ ⊥ i I t + µ k i z z (here we understand these parameters as relativeto the vacuum permittivity ε and the permeability µ , respectively, with I t being the unity dyadic in the transverseplane), the Z -parameters are (see, e.g., Ref. 11): Z i = Z i = j Z TE , TM i sin( β TE , TM i d i ) , Z i = Z i = − j Z TE , TM i tan( β TE , TM i d i ) , (12)where β TE , TM i d i is the electric thickness of the i -th layer, and the wave impedances of a spatial harmonic with wavevector k = ( k x , k y , β TE , TM i ) are Z TE i = η k µ ⊥ i β TE i , Z TM i = η β TM i k ε ⊥ i , (13)where the propagation constants for the two polarizations are expressed through q = q k x + k y as (see, e.g., 11): β TE i = vuut µ ⊥ i ( k ε ⊥ i µ k i − q ) µ k i , β TM i = vuut ε ⊥ i ( k µ ⊥ i ε k i − q ) ε k i . (14)Thus, when considering plane waves of fixed polarization and fixed transverse wave number q ≡ | k t | such a slab isdescribed by a 2 × Z -parameters, much like a four-pole network in the circuit theory. To make theanalogy complete, we may introduce the effective “currents” flowing into this four-pole network and relate them tothe magnetic fields at the two interfaces of the slab as I i, TM1 , = p A k t | k t | · ( n , × H i , t ) , I i, TE1 , = p A z × k t | k t | · ( n , × H i , t ) , (15)and the effective “voltages” at the input and the output interfaces of the slab V i, TM1 , = p A k t | k t | · E i , t , V i, TE1 , = p A z × k t | k t | · E i , t . (16)The factor √ A where A is the unit area in the transverse plane ensures that the complex power V i, TE , TM1 , I i, TE , TM1 , ∗ = − A n , · (cid:16) E i, TE , TM1 , × H i, TE , TM1 , ∗ (cid:17) is trivially related to the complex Poynting vector of a mode. Then, for a givenpolarization (TE or TM) relation (11) assumes the form (cid:18) Z i Z i Z i Z i (cid:19) · (cid:18) I i I i (cid:19) − (cid:18) V i V i (cid:19) = (cid:18) e i e i (cid:19) . (17)The equivalent circuit that corresponds to this relation is shown in Fig. 2. The two EMFs at the input and theoutput of this circuit represent the effect of thermal fluctuations inside the chosen layer. The mean square amplitudeof these equivalent sources is derived in Appendix using the approach of the distributed fluctuating current. However,one can apply the Nyquist theory directly to the electric circuit in Fig. 2 and obtain the same result. Namely, bydisconnecting the load from the output of the four-pole network, i.e., setting I i = 0, we eliminate the contribution of e i and obtain a simple two-pole network with the input impedance Z in = V i /I i = Z i . Therefore, the mean squareamplitude of the fluctuating EMF e i within an angular frequency interval ∆ ω is( e i ) = 2Θ( ω, T i )Re( Z i ) ∆ ωπ (18)Repeating the same procedure while interchanging the roles of the input and the output one obtains that( e i ) = 2Θ( ω, T i )Re( Z i ) ∆ ωπ (19) FIG. 2: Equivalent four-pole network of a material layer under temperature T = T i .FIG. 3: (Color online) An illustration to the calculation of the radiative thermal flux through a selected boundary z = z i in amultilayered structure. It is evident that, in general, the fluctuating sources e i and e i must be partially correlated, because they bothrepresent the fluctuations within the same layer. Therefore, in calculations involving expressions which are quadraticin voltage and (or) current (e.g., power) one may also need the correlation function of these EMFs: ( e i e i ). As isshown in Appendix, this correlation can be presented as:( e i e i ) = 2Θ( ω, T i )Re( Z i ) ∆ ωπ . (20)The set of relations (18)–(20) can be also obtained directly from the FDT. Indeed, if one identifies the charges q i , = I i , / ( jω ) as the state variables of the circuit depicted in Fig. 2 and the EMFs e , as the random forcesassociated with the fluctuations, then the FDT demands that for the fluctuations concentrated within a narrowfrequency interval ∆ ω ( e in e im ) = j ~ (cid:2) ( α − ) ∗ mn − ( α − ) nm (cid:3) coth ~ ω k B T × ∆ ωπ , (21)where α mn are the generalized susceptibilities such that q m = P n α mn ( ω ) e n . It is readily seen that ( α − ) mn = jωZ mn .Substituting this into (21) while taking into account the symmetry properties of Z mn we obtain (18)–(20) afterdropping the irrelevant contribution resulting from the quantum zero-point fluctuations.The relations (18)–(20) together with (17) written for both polarizations fully describe the fluctuations within amaterial layer. A structure formed by many layers can now be equivalently represented by a chain connection of manyfour-pole networks each representing a layer. Let us now select an arbitrary boundary between a pair of layers in amultilayered structure and find the radiative power flux per unit area of this boundary (Fig. 3).Because the fluctuating EMFs belonging to separate layers are uncorrelated, we may first consider only the sourceswhich are located at z < z i , where z i is the position of the selected boundary. We number the layers and theboundaries such that the i -th layer is located at z i − ≤ z ≤ z i . Then, the layers in the half-space z > z i can beconsidered as passive (no radiation is comming out of them). In the circuit theory terms these layers constitute aload for the other, active, part of the structure located at z < z i , and can be equivalently represented by an inputimpedance, which, for a chain connection of four-pole networks, is given by the following recursive formula Z i +1in+ = Z i +111 − Z i +112 Z i +121 Z i +122 + Z i +2in+ , (22)where Z i +2in+ is the input impedance of all the layers behind the ( i + 1)-th layer. The recursion is terminated withthe input impedance of the last layer which extends up to z = + ∞ (it can be, for instance, free space behind thestructure), i.e., with the wave impedance of the last layer. Substituting (12) into (22) for the two main polarizationsin uniaxial layers we obtain Z i +1in+ = Z i +1 Z i +2in+ + jZ i +1 tan( β i +1 d i +1 ) Z i +1 + jZ i +2in+ tan( β i +1 d i +1 ) , (23)where, for brevity, Z i +1 ≡ Z TE , TM i +1 and β i +1 ≡ β TE , TM i +1 . Eq. (23) is well-known in the theory of transmission lines.On the other hand, we may apply Th´evenin’s theorem to the active layers located at z < z i . Doing this, one firstfinds the internal impedance of Th´evenin’s equivalent circuit: Z i in − = Z i − Z i Z i Z i + Z i − − , (24)or, after substituting (12), Z i in − = Z i Z i − − + jZ i tan( β i d i ) Z i + jZ i − − tan( β i d i ) , (25)where Z i − − is the internal impedance of the rest of the layers located at z < z i − , and the recursion terminates atthe layer which extends down to z = −∞ . Note that because here we consider reciprocal structures, the Th´eveninimpedance Z i in − equals the input impedance of all the layers located at z < z i as seen by a wave incident from thehalf-space z > z i .Next, the equivalent voltage generator in Th´evenin’s theorem (recall that the EMF of this generator is the same asthe output voltage of the network under the open circuit condition) can be found recursively as: E i g = − e i + Z i Z i + Z i − − ( e i + E i − ) , (26)where E i − is the equivalent EMF of all the sources located at z < z i − . This EMF is defined at the boundary z = z i − . Taking into account relations (18)–(20) and the fact that the EMFs corresponding to distinct layers are notcorrelated, the mean square amplitude of fluctuations of E i g can be expressed after some algebra as( E i g ) = ( e i ) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z i Z i + Z i − − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( e i ) − Z i Z i + Z i − − ! ( e i e i ) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z i Z i + Z i − − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( E i − ) = 2 R i th Θ( ω, T i ) ∆ ωπ + F i ( E i − ) , (27)where F i = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z i Z i + Z i − − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = | Z i | | Z i cos( β i d i ) + jZ i − − sin( β i d i ) | , (28)and R i th = Re( Z i in − ) − F i Re( Z i − − ) . (29)This important result shows that the effect of thermal fluctuations within the i -th material layer under the temperature T i is fully equivalent to the effect of fluctuations in a resistance R i th placed under the same temperature.In other words, Eq. (29) manifests that the input resistance Re( Z i in − ) of a stack of layers can be split into twoaddends: Re( Z i in − ) = R i th + F i Re( Z i − − ). When considered together, these addends represent the total loss in thestack. However, when the thermal fluctuations are of concern, Eq. (29) allows us to separate explicitely a part of theinput resistance that appears in the Nyquist formula as being under physical temperature of the i -th layer. Thus, thenoise produced within the i -th layer is associated with R i th . The other addend, F i Re( Z i − − ), is due to the loss in thelayers located below the i -th layer, and the thermal noise associated with it is understood as the noise received by the i -th layer from the background. FIG. 4: Th´evenin’s equivalent network of thermal noise in layered media.
Similar concepts exist, for example, in the antenna theory where the thermal noise of an antenna is represented asa sum of the noise generated locally by the ohmic loss in the antenna (analogous to R i th ) and the noise received fromthe environment. The first addend in this case is proportional to the antenna loss resistance (which vanishes for anantenna made of a perfect conductor) and the second term is proportional to the radiation resistance of the antenna.From Eq. (29), Re( Z i − − ) = R i − + F i − Re( Z i − − ), therefore, we may as well write R i th = Re( Z i in − ) − F i R i − − F i F i − R i − − · · · , (30)where the series terminates at the layer (with the index i − M ) that extends to z = −∞ , for which R i − M th =Re (cid:0) Z i − M in − (cid:1) = Re (cid:16) Z TE , TM i − M (cid:17) . We may analogously expand the last addend in (27) which corresponds to the effect offluctuations in the layers located at z < z i − . Doing so we obtain( E i g ) = 2 R i th Θ( ω, T i ) ∆ ωπ + 2 F i R i − Θ( ω, T i − ) ∆ ωπ + 2 F i F i − R i − Θ( ω, T i − ) ∆ ωπ + · · · , (31)Thus, we conclude that the effect of thermal fluctuations in all layers located at z < z i is the same as in a chainof resistors with the values R i eff = R i th , R i − = F i R i − , R i − = F i F i − R i − , etc., kept under the temperatures T i , T i − , T i − , etc. This result is analogous to the known formula for the thermal noise in cascaded amplifiers, in whichcase the quantities F i are called the noise factors.The corresponding equivalent circuit is shown in Fig. 4, in which we split Th´evenin’s internal impedance into areactive part X gen ≡ Im( Z i in − ) and a resistive part Re( Z i in − ) = P n R i − n eff . Respectively, Th´evenin’s EMF splits into aseries of uncorrelated fluctuating EMFs: E i g = P n E i − n eff , with ( E i − n eff ) = (2 /π )Θ( ω, T i − n ) R i − n eff ∆ ω , n = 0 , , , . . . , M, representing the effect of thermal fluctuations in the layers located at z < z i . The rest of the structure at z > z i ismodeled by an effective load impedance Z load ≡ Z i +1in+ .The radiative heat flux from any layer located at z < z i into the half-space z > z i can now be trivially calculatedbased on this equivalent circuit. Namely, the power spectral density associated with the plane waves with a giventransverse wave vector k t = ( k x , k y ) and a given angular frequency ω , originated in the slab with the index i − n , canbe expressed as P i − nω, k t = 1∆ ω ( E i − n eff ) | Z i in − + Z load | Re( Z load ) = 2 π Θ( ω, T i − n ) R i − n eff | Z i in − + Z i +1in+ | Re (cid:0) Z i +1in+ (cid:1) . (32)Respectively, for the total radiative heat flux (associated with waves of a selected polarization: TE or TM) into thehalf-space z > z i we have S z>z i = X n ∞ Z dω Z Z dk x dk y (2 π ) P i − nω, k t = 12 π X n ∞ Z dω ∞ Z qdq P i − nω, k t . (33)The flux crossing the the same boundary in the opposite direction is found by reversing the roles of the active andpassive layers. IV. PARTICULAR CASE I: BLACK BODY RADIATION
Let us apply this equivalent circuit theory to calculate the power radiated by a black body per unit of frequencyand surface area. We assume that a very thick black body lies in the lower half-space z < T . The upper half-space z > z = 0.The equivalent circuit for this system is composed of a single resistance R (1)eff = Re (cid:16) Z (1)in − (cid:17) , where Z (1)in − is the inputimpedance of the half-space z < E (1)eff , a reactance X gen = Im (cid:16) Z (1)in − (cid:17) , and a load Z load = Z (2)in+ = Z TE , TM0 , which is the input impedance of the open half-space at z > Z (1)in − = Z TE , TM0 with Z TE0 = p µ /ε p − q /k , Z TM0 = r µ ε q − q /k , (34)and, using Eq. (32), we may write the power spectral density associated with the radiative heat flux into the openhalf-space as P ω, k t = 2Θ( ω, T ) π h Re (cid:16) Z TE , TM0 (cid:17)i (cid:12)(cid:12)(cid:12) Z TE , TM0 (cid:12)(cid:12)(cid:12) = (cid:26) Θ( ω, T ) / (2 π ) , q ≤ k , q > k (35)The result is the same for both polarizations. From here, the total power emitted from the black body by bothpolarizations per unit of its surface, per unit of frequency is dSdω = 2 × Z Z dk x dk y (2 π ) P ω, k t = Θ( ω, T )4 π Z Z k x + k y In this section we derive a generalization of the Polder-van Hove formula applicable to layered uniaxial magneto-dielectrics, using the equivalent circuit model developed in Section III. The geometry of the structure is the same as inFig. 1 (a). We are interested in the radiative thermal transfer between the media which occupy the half-spaces z < z > d (the media with indices 1 and 3). These half-spaces are kept under temperatures T and T , respectively.The region 0 < z < d (the region 2) may be filled with another uniaxial medium kept under temperature T , or maybe left empty (which is the case of a vacuum gap).In order to find the radiative heat flux from medium 1 into medium 3 we split the structure at the plane z = d , andconsider the layers located at z < d as active layers. The half-space z > d plays the role of a load. The equivalentcircuit of such a structure can be represented as in Fig. 4, with a pair of resistors R (1)eff = F (2) R (1)th = F (2) Re( Z ) and0 R (2)eff = Re (cid:16) Z (2)in − (cid:17) − R (1)eff , a pair of the corresponding fluctuating EMFs E (1)eff and E (2)eff , a reactance X gen = Im (cid:16) Z (2)in − (cid:17) ,and a complex load Z load = Z .One may verify that in the vacuum gap case R (2)eff = 0, which is a consequence of the fact that there is no dissipationin the gap. Evidently, the same conclusion holds when the gap is filled with a lossless medium. Nevertheless,the following derivation is general enough to be applicable to both medium-filled or vacuum gaps, with or withoutdissipation.We start with calculating the noise factor F (2) . Because Z (1)in − = Z , we obtain from (28): F (2) = | Z | | Z cos( β d ) + jZ sin( β d ) | = | − Γ | | e − jβ d | | − Γ e − jβ d | , (40)where Γ is defined in Section II. Thus, from (32), P → ω, k t = 2 π Θ( ω, T ) F (2) Re( Z ) | Z (2)in − + Z | Re( Z ) . (41)Next, Th´evenin’s internal impedance Z (2)in − is found from (25): Z (2)in − = Z Z + jZ tan( β d ) Z + jZ tan( β d ) = Z e − jβ d − Γ e − jβ d , (42)from which the total impedance of the network reads Z (2)in − + Z = 2 Z − Γ − Γ Γ e − jβ d − Γ e − jβ d , (43)where Γ is defined in Section II. Substituting (40) and (43) into (41) we obtain P → ω, k t = Θ( ω, T )2 π | − Γ | | − Γ | | e − jβ d | | − Γ Γ e − jβ d | Re( Z )Re( Z ) | Z | . (44)Respectively, the total radiative heat flux associated either with TE or TM polarized waves originated from medium1 and absorbed in medium 3 is S → = 1 π ∞ Z dω ∞ Z qdq Θ( ω, T ) N ( ω, q ) , (45)where N ( ω, q ) = | − Γ | | − Γ | | e − jβ d | | − Γ Γ e − jβ d | Re( Z )Re( Z )4 | Z | . (46)In the vacuum gap case, the wave impedance Z and the propagation factor β are purely real (imaginary) for thepropagating waves (evanescent waves) in the gap. Therefore, because Z i = Z (1 + Γ i ) / (1 − Γ i ), i = 1 , 3, we haveRe( Z i ) = | Z | − | Γ i | | − Γ i | , (propagating waves), −| Z | i ) | − Γ i | , (evanescent waves), (47)which results in the Polder-van Hove formulas when substituted into (45) and (46). However, the more general resultrepresented by Eqs. (45)–(46) holds for arbitrary uniaxial magneto-dielectric media filling the gap. It is easy to verifythat Eq. (46) can be as well written in form (9).Moreover, in general, when the gap is filled with a lossy medium and T = 0 one also has to take into accountthe radiative heat flux to medium 3 that is originated in the gap (i.e., in medium 2). The power spectral densityassociated with it is, from Eq. (32), P → ω, k t = 2 π Θ( ω, T ) R (2)eff | Z (2)in − + Z | Re( Z ) = Θ( ω, T )2 π | − ˜Γ | | − Γ | || − ˜Γ Γ | R (2)eff Re( Z ) | Z | , (48)1where ˜Γ = Γ e − jβ d is the reflection coefficient defined at the plane z = d . One may note that (48) has thesame form as (44) with d = 0 and one of the wave impedances replaced by the effective resistance R (2)eff . This is theconsequence of the fact that when the thickness of the middle layer increases, ˜Γ → R (2)eff → Re( Z ), and (48)reduces to the Polder-van Hove’s result for two media in direct contact. VI. ANTENNA THEORY AND CIRCUIT THEORY CONCEPTS APPLIED TO RADIATIVE HEATTRANSFER In this section, in order to better understand the circuit model developed in this work, we establish a connectionbetween our model and the classical theory of noise in receiving antennas. We consider the radiative heat transferexample from Section V and employ an analogy between the heat-receiving half-space (medium 3) and a loadedreceiving antenna. More exactly, in this analogy the unit area of the interface between the media 2 and 3 is treated asan aperture antenna that receives the power of thermal radiation from the half-space z < d and delivers it to medium3 which is understood as the antenna load. In what follows, we assume that medium 2 is lossless, therefore, all theradiative heat delivered to medium 3 is generated in medium 1.The equivalent circuit of this problem is the same as the one discussed in Section V with R (2)eff = 0. Thus, in theantenna analogy there is one noise source with the internal impedance Z (2)in − ≡ Z A ≡ R A + jX A , where X A = X gen isthe antenna reactance and R A = F (2) Re( Z ) is the radiation resistance of the antenna. Such analogy signifies thatthe effective radiation resistance of an aperture equals the real part of the input impedance of the half-space seenfrom the aperture. The antenna load is represented in this analogy by the impedance Z load = Z which is the waveimpedance of medium 3.As is known from the theory of noise in lossless antennas (an aperture by itself has no loss), the mean-square EMF e of the thermal noise of a directive antenna (e.g., a radio telescope) is given by the Nyquist formula (1), in whichone inserts the radiation resistance of the antenna [as R ( ν )] and the effective temperature of the area of the sky towhich the antenna is directed (which is medium 1 under temperature T = T in our analogy). Applying this to theantenna equivalent circuit, we may write for the noise power at the antenna load: P out = e | Z A + Z | Re( Z ) = 2 π Θ( ω, T ) R A | Z A + Z | Re( Z ) . (49)For our example of a lossless medium 2, R A = Re (cid:16) Z (2)in − (cid:17) , which brings us to the same result as the more generalcascade-circuit model (41), i.e., P out = P → ω, k t . Let us also note that if there would be no separating layer, then Z A would be simply equal to the wave impedance of medium 1. With the separation layer in place, antenna model calcu-lation is equivalent to calculation of input impedance of a transmission-line section loaded with a known impedance[Eq. (42)].It is worth noting that the tight connection between the model of the present paper and the theories of noise inantennas and cascaded electric networks allows one to better understand optimal conditions for radiative heat transferthrough composite layers and, consequently, to design these material structures aiming for desired and optimizedperformance. In particular, from (49) and (41) we see that the problem of maximizing radiative heat transfer forgiven layer temperatures reduces to an equivalent problem of matching a generator to a load. Let us discuss this issueassuming for simplicity, that the media 1 and 3 are the same, i.e., Z = Z . At the first glance, it appears that theoptimal heat transfer is ensured if we simply connect the two equivalent media together (or fill the gap with the samemedium as those 1 and 3). However, this is true only if the wave impedance is real. For complex Z and Z (whichis realistic even for propagating modes in view of losses in the media), the best radiative heat transfer correspondsto the conjugate impedance matching Z = Z ∗ when the negative reactance of one of two media is compensated bythe positive reactance of the other one. Then, the spatial spectrum of the heat transfer function in accordance to (9)turns to N = 1 / Z = R + jX we have N = R / | Z | . This tells us that the radiative heat transfer through a properly filled gap can be in principle madelarger than that through the direct contact of two equivalent media.In Ref. 5 it was proposed to insert a slab of the so-called indefinite medium between media 1 and 3 [Fig. 1 (b)].Indefinite media (also called hyperbolic metamaterials ) are uniaxial dielectrics characterized by the permittivitytensor that has opposite signs of the longitudinal and transverse components. Such filling allows for enhancementof the heat transfer by increasing the noise factor F (2) . This factor increases compared to the vacuum gap becauseindefinite media support propagation of spatial harmonics with high transverse wavenumbers, which would otherwisebe evanescent in the gap [notice the exponential factor in Eq. (40)]. For such waves, F (2) dramatically increases at thespatial frequencies that correspond to the minima of the denominator of Eq. (40). However, it is not straightforward2to ensure proper impedance match in such structures. In this view, the structure suggested and studied in Ref. 5is not fully optimal, although it still demonstrates that a specifically crafted filling may dramatically enhance theradiative heat transfer through the gap.The transmission line analogy, however, suggest an immediate possibility to circumvent the problem of impedanceconjugate match. The key is to make the wave impedances of all three media real (at least, approximately) andequal in all layers. This can be realized using a wire medium in region 2 which extends inside regions 1 and 3.Indeed, it is known that wire media support propagating TEM modes with high spatial frequencies (transverse wavenumbers) q , including q > k and limited only by the period of the wire array (see, e.g., in Ref. 23). Wires should begood conductors at the frequency of interest, and the background media nonconducting. If the wires extend over allthree regions and none of them has a conducting background, then these TEM modes will exist and have real waveimpedance everywhere in the system, in principle allowing the conjugate match of the load to the heat source for awide range of spatial frequencies q .In the present paper we, however, do not present the studies of the heat transfer optimization keeping these resultsfor our next publications. Instead, in the next section we report some numerical results illustrating the applicabilityof our model for solving practical problems. VII. NUMERICAL EXAMPLE: RADIATIVE THERMAL TRANSFER THROUGH ANANOSTRUCTURED LAYER In order to demonstrate the applicability of the developed theory to real-world problems, we calculate in this sectionthe spectral density of the radiative heat flux absorbed in medium 3 of the structure depicted in Fig. 1 (b): s ≡ dS → dω = 12 π ∞ Z P → ω, k t qdq. (50)We compare the heat transferred across the vacuum gap with that transferred across the gap filled with eithernormally oriented metal-state single-wall carbon nanotubes (CNT) or with similarly oriented golden (Au) nanowires.Calculation for the array of CNT is done in order to validate the present model using the exact simulations of Ref. 5.Calculations for metal nanowires are done in order to confirm or decline the effect of giant enhancement of radiativeheat transfer in the near infrared (IR) range due to the presence of nanowires. This effect was predicted but notstudied in Ref. 5. In accordance with the theory presented above, all calculations are done for homogenized media.The homogenization models for i -th medium results in explicit formulas for ε ⊥ i and ε k i (whereas µ ⊥ i = µ k i = 1).The mid-IR homogenization model for the array of CNT was described and validated in Ref. 24. The parametersof the array of CNT correspond to those of Ref. 5 [see also in Fig. 1 (b)]. The homogenization model for alignedmetal nanowires representing their array as a layer of an indefinite material was described in Ref. 25. It is applicableto the visible and near-IR ranges where it offers a rather high accuracy for optically dense arrays of rather thin wires.Practically, for the band λ =1–2 µ m one needs the period below 300–600 nm and the wire diameter of the order ofthe skin depth in the metal or smaller. For Au in this wavelength range it means that the thickness of the wiresshould not exceed 20–50 nm.In both cases (CNT in the mid-IR range and nanowires in the near IR) the presence of the metamaterial enhancesthe radiative heat transfer. This effect results from the conversion of TM-polarized evanescent waves into propagatingones in the effective indefinite material filling the gap. Respectively, in this section we analyze only that part ofthe radiative heat which is transferred by TM-polarized waves.In this numerical example we neglect the contribution of thermal sources located in medium 2 (i.e., in CNT andnanowires). The impact of thermal production and absorption in medium 2 will be evaluated in our next paper, wherewe will also consider thermo-photovoltaic applications of our present model. As is clear from Fig. 1 (b), medium 2 inthe gap is a triple-layered structure, therefore, there are in total five material layers in the whole structure. However,since in this example we neglect the thermal processes in the gap, we may replace the layers in the gap by a singleequivalent four-pole network. Its Z -matrix, Z (2) mn , is obtained in a standard manner from its transfer matrix, withthe latter being a product of transfer matrices of effectively homogeneous anisotropic layers with thicknesses h , h and h . Next, the power spectral density P → ω, k t in Eq. (50) is calculated from Eq. (41) where the noise factor is F (2) = (cid:12)(cid:12)(cid:12) Z (2)21 (cid:12)(cid:12)(cid:12) / (cid:12)(cid:12)(cid:12) Z (2)11 + Z (cid:12)(cid:12)(cid:12) [Eq. (28)], and the internal impedance is Z (2)in − = Z (2)22 − Z (2)12 Z (2)21 / (cid:16) Z (2)11 + Z (cid:17) [Eq. (24)].In Fig. 5 (a) we depict the energy transfer coefficient | τ | introduced in Ref. 5 [Eq. (12) of Ref. 5] which differsfrom our heat transfer spatial spectrum N ( ω, q ) defined by Eq. (10), by the factor 4 | Z | /R R . In the present casemedia 1 and 3 are equivalent (heavily doped silicon), i.e., Z = Z , R = R . All parameters in the calculation3 π | τ | CNT, modelVacuum gap, modelVacuum gap, exactCNT, exact (a) N ( q , λ = . µ m ) Nanowires in media 1 and 2Nanowires in medium 2Nanowires in medium 1No nanowires (b) FIG. 5: (Color online) (a) – Coefficient | τ | = 4 | Z /R | N at λ = 7 . µ m versus normalized spatial frequency qa/π for the gapfilled with free space and for that filled with an array of CNT. Dashed lines – exact (beyond homogenization) simulations fromRef. 5. Solid lines – analytical calculations in accordance to the present model. (b) – Spatial spectrum N of the radiative heattransfer function versus dimensionless spatial frequency q/k at λ = 1 . µ m. Calculations are done for four cases – nanowiresare only in medium 2, nanowires are in both media 1 and 2, nanowires are in medium 1 only and nanowires are absent. G a i n , d B Nanowires in media 1 and 2Nanowires in medium 2 (a) G a i n , d B Nanowires in media 1 and 2Nanowires in medium 2 (b) FIG. 6: (Color online) The gain in the spectral density of the transferred radiative heat flux due to the presence of nanowires(two cases of their arrangement). (a) – Gap of thickness d = 2 µ m. (b) – Gap of thickness d = 0 . µ m. illustrated by Fig. 5 (a) correspond to those from Ref. 5. The value | τ | is calculated at wavelength λ = 7 . µ mas a function of the normalized spatial frequency qa/π , where a = 20 nm is the period of the CNT array in thedomains h = d/ h = d/ 3, where the gap thickness is d = 1 µ m. In the domain h = d/ a = √ 20 nm. Heavily doped silicon supports so-called surface plasmon-polariton (SPP) waves generatedon the surfaces of the Si half-spaces at ( qa/π ) = 0 . 01 in the case when the gap is empty. The value q = 0 . π/a nearly corresponds to q = 1 . k , where k is the free-space wave number. The manifestation of this SPP is themaximum of the corresponding curve in Fig. 5 (a). Dashed curves in Fig. 5 (a) correspond to exact simulations ofRef. 5 which take into account the microstructure of the material layer in the gap. In Fig. 5 (a) we show only theregion of spatial frequencies ( qa/π ) ≤ . 08 in which the difference between the exact and homogenized models of theCNT array is negligibly small (see Ref. 5). For both empty gap and gap filled with CNT the agreement betweenour model and the exact simulations is excellent. Local maxima of the transmittance spatial spectrum for the gap4filled with CNT correspond to the thickness resonances of spatial harmonics [in presence of CNT the whole region( qa/π ) ≤ . 08 corresponds to propagating plane waves, though the inequality q > k holds for ( qa/π ) > . · − ].Fine agreement between our circuit model and exact simulations pertains at other wavelengths besides the onesmentioned in Figs. 5 (a,b), because with our circuit model we also have reproduced the frequency dependence of thegain in the heat transfer G ( ω ) = s ( CNT )13 ( ω ) /s (0)13 ( ω ), calculated in Ref. 5. Here s (0)13 corresponds to the vacuum gapand s ( CNT )13 corresponds to the gap filled with CNT.In Fig. 5 (b) we present the dependence N ( q, λ ∗ ) for the case when the gap d = 2 µ m is filled with golden nanowires.The complex permittivity of gold in the range λ = 1 . . . µ m was taken from Ref. 26. Regions 1 and 3 in this caseare filled with doped germanium used in real thermo-photovoltaic systems whose complex permittivity was takenfrom Ref. 27. Function N ( q, λ ∗ ) was calculated for nanowires with volume fraction p = 0 . h and h and p = 0 . h (when h = h = h ). Here λ ∗ = 1 . µ m has been chosen having in mind possiblethermophotovoltaic applications (at this wavelength the doped Ge has nearly maximal photovoltaic spectral response).The result for N (thick solid curve) was compared with that for the empty gap d = 2 µ m (thin dashed curve) at thesame wavelength.It has to be mentioned that the structure shown in Fig. 1 (b) with free-standing metal nanowires is an abstraction— in a feasible structure nanowires are partially submerged into the host material. Therefore we have also calculatedthe function N ( q, λ ∗ ) for the case when Au nanowires are semi-infinite and have the same volume fraction p = 0 . h = h = 0 and h = d . In this case the nanowres touch the surface of medium 3 andthe thermal transfer by the direct thermal conductance may be of significance, in addition to the radiative one. Thiseffect is not considered in the present paper. Additionally, we have studied the case when nanowires with p = 0 . N ( q ) with respect to the empty gap is very significant forboth cases when the nanowires fill in the gap. Function N ( q ) in the case of the vacuum gap has two local maxima inthe region q < k resulting from Fabry-Perot resonances. Because the structure is not fully impedance-matched thesemaxima are much smaller than the achievable limit N = 1 / N ( q ) vanishes fast at q > k . Unlike the situationillustrated by Fig. 5 (a) the real part of the complex permittivity of Ge at λ = 1 . µ m is positive and SPP cannotbe excited. Since the value of N ( q ) for the empty gap is so small, the gain granted by Au nanowires to the heattransfer between two half-spaces of Ge turns out to be larger than that offered by CNT to the heat transfer betweentwo half-spaces of Si calculated in Ref. 5.The presence of nanowires only in medium 1 turns out to be destructive for the amplitude of N . In this casemedium 1 is an indefinite metamaterial, and the mismatch between medium 1 and free space increases. However, ifnanowires are present in both media 1 and 2, noticeable values of N keep for q > k . The case when nanowires arelocated only in medium 2 is the best one: it corresponds to the smallest mismatch between the media.The dependencies shown in Fig. 5 (b) are typical for every λ in the band of the photovoltaic operation of Ge( λ = 1 − µ m). As a result, due to the presence of nanowires the heat transfer gain G = s ( NW )13 /s (0)13 is almost uniformover a wide range of wavelengths. Here s ( NW )13 corresponds to the case of nanowires in the gap and s (0)13 corresponds tothe vacuum gap. In Fig. 6 (a) we present the gain in dB, i.e., 10 log G calculated for the case d = 2 µ m. The hugegain keeps for the interval of values d = 0 . . . . µ m. In Fig. 6 (b) we show the same gains for the case d = 0 . µ m.We can conclude that the presence of Au nanowires in the micrometer or submicron gap can increase the near-IRenergy transfer across the gap by 3 orders of magnitude. This result confirms the expectations of Ref. 5. VIII. CONCLUSIONS In this work we have formulated an equivalent circuit theory of the radiative heat transfer in uniaxial stratifiedmagneto-dielectric media. We have proven that the effect of thermal-electromagnetic fluctuations in such structurescan be fully determined without an explicit knowledge of the microstructure of the layers, as well as without a needto employ any calculations based on distributed fluctuating currents. Instead, the only physical characteristic onwhich we base our theory is the effective input impedance of a stack of layers, which can be obtained for any spatialharmonic of the field (including both propagating and evanescent waves) using the methods of VTLT. We have shownthat such impedance representation, while being in full agreement with sophisticated full-wave methods known fromthe literature, results in simple formulas analogous to Nyquist theory-based formulas for thermal noise in cascadedelectric circuits (for example, cascaded amplifiers). Therefore, with this model some important concepts from thetheory of electric networks (conjugate-impedance match, optimal filtering, etc.) can be imported into the field ofradiative thermal transfer in multilayered structures.5From the point of view of practical implementations, the developed equivalent circuit approach offers significantsimplifications as compared to the known theories of radiative heat transfer based on distributed fluctuating currents.Without any modifications, our method can be used in heat transfer studies in uniaxial anisotropic media that includemicron and (or) submicron-thick layers. Moreover, our model is applicable to radiative heat transfer in compositeor nanostructured layers (if these layers are effectively homogeneous for spatial harmonics of the electromagneticfield which transfer the radiative heat), and is readily generalizable to stratified bi-anisotropic and spatially dispersivematerials. Therefore, we hope that our work may significantly enlarge the scope of the radiative heat transfer researchin composites, especially in nanostructured metamaterials. We believe that this may lead to new opportunities in thedesign of efficient thermal energy harvesting devices, like thermophotovoltaic converters and such. Appendix: Nyquist formula for a reciprocal anisotropic and lossy magneto-dielectric slab We consider a uniaxial magneto-dielectric slab described by the macroscopic Maxwell equations for the time-harmonic fields ∇ × E = − jωµ a · H − J m , ∇ × H = jωε a · E + J e , (51)with the absolute permittivity and permeability dyadics of the form ε a = ε ( ε ⊥ I t + ε k z z ) and µ a = µ ( µ ⊥ I t + µ k z z ). We assume that the material of the slab is lossy, therefore, by the fluctuation-dissipation theorem (atnon-zero temperature) there appear fluctuating external currents J e and J m in the slab. The explicit form of thesecurrents is not important at this stage. As in the main text, here we use the convention in which the time-harmonicquantities are understood as root mean square (rms) values.Let the fields E ′ , H ′ be an arbitrary solution of the Maxwell equations (51) with J e = J m = 0 within the slab. Then,considering the two systems of Maxwell equations with non-zero sources and with vanishing sources, respectively, wecan form the Lorentz lemma ∇ · ( E × H ′ − E ′ × H ) = E ′ · J e − H ′ · J m . (52)Integrating it over the volume V of the slab, we obtain the reciprocity relation Z S n · ( E × H ′ − E ′ × H ) (cid:12)(cid:12) S dS + Z S n · ( E × H ′ − E ′ × H ) (cid:12)(cid:12) S dS = Z V ( E ′ · J e − H ′ · J m ) dV, (53)where S , are at the two interfaces of the slab, and n , are the outer unit normals to these surfaces, respectively.One may select any solution of the uniform Maxwell equations within the slab for the fields E ′ , H ′ . For us it isconvenient to use the one that has the form E ′ ( r ) = E ′− k t ( z ) e j k t · r , H ′ ( r ) = H ′− k t ( z ) e j k t · r , (54)where the z -axis is orthogonal to the slab and the real vector k t lies in the plane of the slab (the xy plane). Physically,such a form corresponds to a superposition of plane waves with the same transverse wavenumber: − k t . In Eq. (54), E ′− k t and H ′− k t define the field solution profile within the slab as a function of z .Substituting (54) into the reciprocity relation (53), we obtain (the first slab interface is at z = z and the secondone is at z = z ) n · ( E k t × H ′− k t − E ′− k t × H k t ) (cid:12)(cid:12) z = z + n · ( E k t × H ′− k t − E ′− k t × H k t ) (cid:12)(cid:12) z = z = Z z z ( E ′− k t · J e k t − H ′− k t · J m k t ) dz, (55)where we have decomposed the fluctuating currents J e,m and the fields E , H into plane waves using the Fouriertransform defined as F ( r ) = A (2 π ) Z Z F k t ( z ) e − j k t · r d k t , F k t ( z ) = 1 A Z Z F ( r ) e j k t · r d r , (56)where A is the unit area in the xy -plane, and F can be any of the fields or currents.Eq. (55) is the reciprocity relation for the wave components characterized with a fixed transverse wavenumber. Inorder to simplify further writing we will use the notation F , ≡ F ± k t ( z , ) with F being any of the fields or currents.6Then, noticing that only transverse components of the fields play any role on the left-hand side of (55) we rewrite itas n · ( E t × H ′ t − E ′ t × H t ) + n · ( E t × H ′ t − E ′ t × H t ) = E ′ t · ( n × H t ) − E t · ( n × H ′ t ) + E ′ t · ( n × H t ) − E t · ( n × H ′ t ) . (57)Let us remind that the quantities E ′ , t and H ′ , t have the meaning of the transverse components of the electricand magnetic fields at the interfaces of a source-free magneto-dielectric slab. Therefore, as follows from the vectortransmission line theory (VTLT) for such slabs, these components are related by the impedance matrix of the slab (cid:18) E ′ t E ′ t (cid:19) = Z Z Z Z ! · (cid:18) n × H ′ t n × H ′ t (cid:19) . (58)The components of this matrix are dyadics that are even in k t : Z mn ( − k t ) = Z mn ( k t ). Also, due to the symmetryand the reciprocity, Z = Z , Z = Z , and Z Tmn = Z mn . Using (58) on the left-hand side of Eq. (57) we obtain A h E ′ t · ( n × H t ) − E t · ( n × H ′ t ) + E ′ t · ( n × H t ) − E t · ( n × H ′ t ) i = I ′ · h Z · I + Z · I − V i + I ′ · h Z · I + Z · I − V i , (59)where we have introduced the vector currents I ′ , ≡ √ A n , × H ′ , t , I , ≡ √ A n , × H , t and the vector voltages V , ≡ √ A E , t , and used the symmetry properties of the impedance dyadics.Let us now work on the right-hand side of (55). At a fixed k t the Maxwell equations for the fields E ′ ( r ), H ′ ( r )reduce to a system of first-order linear differential equations for the vector functions E ′− k t ( z ) and H ′− k t ( z ). From theuniqueness theorem it follows that these functions are univocally defined by boundary conditions imposed either ontangential electric or tangential magnetic field. Thus, we may consider two auxiliary boundary-value problems, thefirst one with the boundary conditions n × H I − k t ( z ) = I ′ / p A , n × H I − k t ( z ) = 0 , (60)and the second one with n × H I − k t ( z ) = 0 , n × H I − k t ( z ) = I ′ / p A . (61)The field equations are the same in these two problems. From linearity it follows that the superposition of the solutionsof the two problems is the same as the fields E ′− k t ( z ) and H ′− k t ( z ) that appear in (59). On the other hand, theseproblems physically correspond to the two cases of the magneto-dielectric slab backed with a magnetic wall (perfectmagnetic conductor, PMC) at z = z and at z = z , respectively.Let us consider the problem with the boundary conditions (60). We may split the vector I ′ into the componentsparallel and orthogonal to k t : I ′ = I ′ , TM k t | k t | + I ′ , TE k t × n | k t | . (62)Thus, the component I ′ , TM corresponds to TM-polarized field, and the component I ′ , TE corresponds to TE-polarizedfield. Because the wave equations in the slab also split into independent equations for the TM and TE waves, we mayalso write for the vector fields E I − k t ( z ) = E ′ , TM ( z ) + E ′ , TE ( z ) , H I − k t ( z ) = H ′ , TM ( z ) + H ′ , TE ( z ) , (63)where the addends are the TM and TE solutions for the fields in the PMC-backed slab. The magnitudes of thesesolutions are proportional to I ′ , TM and I ′ , TE , respectively.Based on the above discussion we find for the right-hand side of (55) Z z z ( E I − k t · J e k t − H I − k t · J m k t ) dz = Z z z ( E ′ , TM · J e k t − H ′ , TM · J m k t ) dz + Z z z ( E ′ , TE · J e k t − H ′ , TE · J m k t ) dz =1 A I ′ · (cid:18) e , TM k t | k t | + e , TE k t × n | k t | (cid:19) , (64)7where e , TM = A I ′ , TM Z z z ( E ′ , TM · J e k t − H ′ , TM · J m k t ) dz, (65) e , TE = A I ′ , TE Z z z ( E ′ , TE · J e k t − H ′ , TE · J m k t ) dz. (66)In an analogous manner we consider the second case with a PMC at z = z and find Z z z ( E II − k t · J e k t − H II − k t · J m k t ) dz = 1 A I ′ · (cid:18) e , TM k t | k t | + e , TE k t × n | k t | (cid:19) , (67)where e , TM = A I ′ , TM Z z z ( E ′ , TM · J e k t − H ′ , TM · J m k t ) dz, (68) e , TE = A I ′ , TE Z z z ( E ′ , TE · J e k t − H ′ , TE · J m k t ) dz. (69)Therefore, combining these results together and using (55), (57), and (59) we obtain I ′ · h Z · I + Z · I − V i + I ′ · h Z · I + Z · I − V i = I ′ · (cid:18) e , TM k t | k t | + e , TE k t × n | k t | (cid:19) + I ′ · (cid:18) e , TM k t | k t | + e , TE k t × n | k t | (cid:19) . (70)Finally, because the vectors I ′ and I ′ are arbitrary, Z Z Z Z ! · (cid:18) I I (cid:19) − (cid:18) V V (cid:19) = (cid:18) e e (cid:19) , (71)where e , = e , , TM ( k t / | k t | ) + e , , TE ( k t × n , ) / | k t | . These equations represent the equivalent vector circuit modelof a magneto-dielectric slab with fluctuating sources. In this model, I , have the meaning of equivalent vector currentsat the two ports of a linear four-pole network of dyadic impedances, and e , are the equivalent vector EMFs actingat the two ports.Because the equivalent EMFs are expressed through the fluctuating currents, they are also fluctuating, stochasticquantities. As is readily seen from (65)–(66) and (68)–(69), the stochastic mean value of the fluctuating EMFs iszero: e , = 0, because J e,m k t = 0. However, the mean-square values of the fluctuating EMFs, as well as their mutualcorrelations are in general different from zero and can be calculated as follows:( e ∗ α,p e β,q ) = A I ′∗ α,p I ′ β,q Z z z ( E ′∗ α,p · J e k t ∗ − H ′∗ α,p · J m k t ∗ ) dz Z z z ( E ′ β,q · J e k t − H ′ β,q · J m k t ) dz ′ = A I ′∗ α,p I ′ β,q Z z z Z z z ( E ′∗ α,p · J e k t ∗ − H ′∗ α,p · J m k t ∗ ) (cid:12)(cid:12) z ( E ′ β,q · J e k t − H ′ β,q · J m k t ) (cid:12)(cid:12) z ′ dz dz ′ = A I ′∗ α,p I ′ β,q (cid:20)Z z z Z z z E ′∗ α,p ( z ) · J e ∗ k t ( z ) J e k t ( z ′ ) · E ′ β,q ( z ′ ) dz dz ′ + Z z z Z z z H ′∗ α,p ( z ) · J m ∗ k t ( z ) J m k t ( z ′ ) · H ′ β,q ( z ′ ) dz dz ′ (cid:21) , (72)where α, β = 1 , p, q = TE , TM. There are no cross terms in the last integral of Eq. (72) because the electricand magnetic fluctuations are statistically independent: the dyadic J e ∗ k t J m k t is such that J e ∗ k t J m k t = 0.8From the fluctuation-dissipation theorem, for the fluctuating currents composed of harmonics within a narrowinterval around a given frequency ω , J e ∗ k t ( z ) J e k t ( z ′ ) = 1 πA jω (cid:16) ε a − ε † a (cid:17) δ ( z − z ′ )Θ( ω, T )∆ ω, (73) J m ∗ k t ( z ) J m k t ( z ′ ) = 1 πA jω (cid:16) µ a − µ † a (cid:17) δ ( z − z ′ )Θ( ω, T )∆ ω. (74)The dimensionality factor 1 /A appears in (73)–(74) because of the form of transformation (56). Note also that (73)–(74) are written for the rms amplitudes of the fluctuating currents.Substituting (73)–(74) into (72) and evaluating the integrals over z ′ we find that( e ∗ α,p e β,q ) = jωA Θ( ω, T )∆ ωπI ′∗ α,p I ′ β,q Z z z h E ′∗ α,p · (cid:16) ε a − ε † a (cid:17) · E ′ β,q + H ′∗ α,p · (cid:16) µ a − µ † a (cid:17) · H ′ β,q i dz. (75)However, from the well-known differential lemma ∇ · ( E × H ∗ + E ∗ × H ) = − jω h E ∗ · (cid:16) ε a − ε † a (cid:17) · E + H ∗ · (cid:16) µ a − µ † a (cid:17) · H i , (76)which holds for arbitrary source-free electromagnetic fields E , ( r ), H , ( r ) within the slab, it follows that Z z z h E ′∗ α,p · (cid:16) ε a − ε † a (cid:17) · E ′ β,q + H ′∗ α,p · (cid:16) µ a − µ † a (cid:17) · H ′ β,q i dz = − jω h n · ( E ′ β,q × H ′∗ α,p ) (cid:12)(cid:12) z = z + n · ( E ′∗ α,p × H ′ β,q ) (cid:12)(cid:12) z = z + n · ( E ′ β,q × H ′∗ α,p ) (cid:12)(cid:12) z = z + n · ( E ′∗ α,p × H ′ β,q ) (cid:12)(cid:12) z = z i , (77)from which we see that if p = q , the integral (77) vanishes due to the orthogonality of the TE and TM polarizations.Next, when p = q and α = β = 1 we obtain from (77), (60)–(61), and (58): Z z z h E ′∗ ,p · (cid:16) ε a − ε † a (cid:17) · E ′ ,p + H ′∗ ,p · (cid:16) µ a − µ † a (cid:17) · H ′ ,p i dz = 2 | I ′ ,p | jωA Re( Z p ) . (78)An analogous result is obtained for α = β = 2. On the other hand, when α = 1 and β = 2 we obtain Z z z h E ′∗ ,p · (cid:16) ε a − ε † a (cid:17) · E ′ ,p + H ′∗ ,p · (cid:16) µ a − µ † a (cid:17) · H ′ ,p i dz = I ′∗ ,p I ′ ,p jωA (cid:0) Z p + Z p ∗ (cid:1) . (79)Combining all these results together and using the reciprocity property of the Z -parameters we find from (75) that( e ∗ α,p e β,p ) = 2 π Re (cid:16) Z pαβ (cid:17) Θ( ω, T )∆ ω, (80)where α, β = 1 , p = TE , TM. Eq. (80) is the generalized Nyquist formula for the thermal-electromagnetic noisein a uniaxial magneto-dielectric layer. H. B. Callen and T. A. Welton, Phys. Rev. , 34–40 (1951). D. Polder and M. van Hove, Phys. Rev. B , 3303 (1971). C. Fu and Z. M. Zhang, Frontiers of Energy and Power Engineering in China , 11 (2007). C. J. Fu and Z. M. Zhang, Front. Energy Power Eng. China , 11 (2009). I. S. Nefedov and C. R. Simovski, Phys. Rev. B , 195459 (2011). S. M. Rytov, Theory of electric fluctuations and thermal radiation , Electronics Research Directorate, Air Force CambridgeResearch Center, Air Research and Development Command, U.S. Air Force, 1959. J. B. Johnson, Phys. Rev. , 87 (1928). H. Nyquist, Phys. Rev. , 110 (1928). D. Pozar, Microwave and RF Design of Wireless Systems , J. Wiley and Sons: NY, 2001, p. 127. T. Kraus, Antennas , McGraw-Hill: NY, 1988, pp. 774-790. S. A. Tretyakov, Analytical modelling in applied electromagnetics , Artech House: Boston-London-Dordrecht, 2003. R. R. A. Syms and L. Solymar, J. Appl. Phys. , 124909 (2011); R. R. A. Syms, O. Sydoruk, and L. Solymar, Phys.Rev. B , 235150 (2011); R. R. A. Syms, L. Solymar, and O. Sydoruk, Proc. Metamaterials’2012, St. Petersburg, Russia,508–510 (2012). F. N. H. Robinson, Noise in electrical circuits , Oxford University: London, 1962. A. van der Ziel, Noise: sources, characterization, measurements , Prentice-hall: Upper Saddle River, NJ, USA, 1970. M. Buckingham, Noise in electronic devices and systems, J. Wiley and Sons, New York, 1983. J. B. Pendry, J. Phys. Cond. Mat. , 6621 (1999). R. Siegel and J. Howell, Thermal radiation heat transfer , 4th ed., Taylor and Francis, New York – London, 2002, p. 525. Zh. Zhang, Nano/microscale heat transfer , McGraw-Hill, Atlanta, Georgia, USA, 2007. D. R. Smith and D. Schurig, Phys. Rev. Lett. , 077405 (2003). E. E. Narimanov and V. Shalaev, Nature , 266 (2007) E. Narimanov, Laser Science, OSA Technical Digest, LWA3 (2011). S.-A. Biehs, M. Tschikin, and P. Ben-Abdallah, PRL , 104301 (2012). C. R. Simovski, P. A. Belov, A. V. Atraschenko, and Yu. S. Kivshar, Advanced Materials , 4229, 2012. I. S. Nefedov, Phys. Rev. B , 155423 (2010). J. Elser, R. Wangberg, E. Narimanov, and V. A. Podolskiy, Appl. Phys. Lett. (2006) 261102. S. Mattei, P. Masclet, and P. Herve, Infrared Physics, , 991 (1989). B. Bitnar, Semiconductor Science Technology , S221 (2003). Throughout this work we use root-mean square (rms) amplitudes for time-harmonic quantities. Thus, if U ( ω ) is the rmsamplitude of u ( t ), then u ( t ) = | U ( ω ) |2