Evaluation of blackbody radiation emitted by arbitrarily shaped bodies using the source model technique

Planck’s famous blackbody radiation law was derived under the assumption that the dimensions of the radiating body are significantly larger than the radiated wavelengths. What is unique about Planck's formula is the fact that it is independent of the exact loss mechanism and the geometry. Therefore, for a long period of time, it was regarded as a fundamental property of all materials. Deviations from its predictions were attributed to imperfections and referred to as the emissivity of the specific body, a quantity which was always assumed to be smaller than unity. Recent studies showed that the emission spectrum is affected by the geometry of the body and in fact, in a limited frequency range, the emitted spectrum may exceed Planck's prediction provided the typical size of the body is of the same order of magnitude as the emitted wavelength. For the investigation of the blackbody radiation from an arbitrarily shaped body, we developed a code which incorporates the fluctuation-dissipation theorem (FDT) and the source model technique (SMT). The former determines the correlation between the quasi-microscopic current densities in the body and the latter is used to solve the electromagnetic problem numerically. In this study we present the essence of combining the two concepts. We verify the validity of our code by comparing its results obtained for the case of a sphere against analytic results and discuss how the accuracy of the solution is assessed in the general case. Finally, we illustrate several configurations in which the emitted spectrum exceeds Planck's prediction as well as cases in which the geometrical resonances of the body are revealed. © 2017 Optical Society of America OCIS codes: (290.6815) Thermal emission; (350.5610) Radiation; (050.6624) Subwavelength structures; (000.4430) Numerical approximation and analysis. References and links 1. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley & Sons, 1983) 2. A. Reiser and L. Schächter, “Geometric effects on blackbody radiation,” Phys. Rev. A 87(3), 033801 (2013). 3. R. Courant and D. Hilbert, Methods of Mathematical Physics (Wiley, 1989), Chap. VI, Sec. 4. 4. H. Weyl, “Das asymptotische Verteilungsgesetz der Eigenwerte linearer partieller Differentialgleichungen (mit einer Anwendung auf die Theorie der Hohlraumstrahlung),” Math. Ann. 71(4), 441–479 (1912). 5. T. Carleman, in Proceedings of the Eighth Scandinavian Mathematics Congress, Stockholm (Ohlsson, Lund, 1935), p. 34. 6. A. Pleijel, “On the eigenvalues and eigenfunctions of elastic plates,” Commun. Pure Appl. Math. 3(1), 1–10 (1950). 7. F. H. Brownell, “An Extension of Weyl’s Asymptotic Law for Eigenvalues,” Pac. J. Math. 5(4), 483–499 (1955). 8. H. P. Baltes and E. R. Hilf, Spectra of Finite Systems (Bibliographisches Institut, 1976). 9. H. B. Callen and T. A. Welton, “Irreversibility and generalized noise,” Phys. Rev. 83(1), 34–40 (1951). 10. L. D. Landau and E. M. Lifshits, Statistical Physics (Pergamon, 1958). 11. S. M. Rytov, Theory of Electric Fluctuations and Thermal Radiation, Bedford, Massachusetts: Air-Force Cambridge Reseach Center, 1959. 12. S. M. Rytov, Y. A. Kravtsov, and V. I. Tatarskii, Principles of Statistical Radiophysics, Vol 3: Elements of Random Fields (Springer-Verlag, 1989). 13. D. Polder and M. Van Hove, “Theory of radiative heat transfer between closely spaced bodies,” Phys. Rev. B 4(10), 3303–3314 (1971). 14. J. J. Loomis and H. J. Maris, “Theory of heat transfer by evanescent electromagnetic waves,” Phys. Rev. B Condens. Matter 50(24), 18517–18524 (1994). Vol. 25, No. 12 | 12 Jun 2017 | OPTICS EXPRESS A589 #285018 https://doi.org/10.1364/OE.25.00A589 Journal © 2017 Received 23 Jan 2017; revised 1 May 2017; accepted 18 May 2017; published 9 Jun 2017


Introduction
Since the dawn of the age of quantum mechanics, Planck's law of blackbody radiation was considered to be one of the cornerstones of modern physics and has usually served as the golden standard for solving problems of thermal radiation. Deviations from this law were attributed to imperfections collectively regarded as emissivity which was always considered to be smaller than unity. However the perception of Planck's law as a basic law of nature has been repeatedly scrutinized in articles and even in textbooks [1]. One of the main limitations of Planck's law stems from the explicit assumption that the radiating cavity is of cuboidal shape whose dimensions are significantly larger than that of the radiated wavelengths. A thorough review of this assumption and examples where it breaks and corrections to Planck's law were presented in [2] and some explicit examples can be found in [3][4][5][6][7][8]. Furthermore, in [2] the main focus is on closed structure, it was emphasized that for proper evaluation of the radiation from an open body, Fluctuation-Dissipation Theorem (FDT) [9] becomes a strict necessity. For the sake of accuracy, we need to point out that in [9] FDT was developed for lumped elements and the relevance to blackbody radiation was assessed based on the electrical dipole impedance. In our study we consider the extension of FDT for distributed systems developed by Landau and Lifshits [10] and Rytov [11][12].
Another assertion that needs to be clarified refers to near vs. far-field contribution. Planck's blackbody radiation formula is relevant only to the far-field contribution and makes no statements regarding the energy stored in near-field modes. Yet it is well known that heat transfer between two macroscopic bodies is dominated by near-field (evanescent) modes if the distance between two macroscopic bodies is of the order of the dominant wavelength of the blackbody radiation [13][14][15][16][17][18][19][20]. Consequently, the conclusions drawn in [1,2] regarding exceeding the Planck's law limit is reasonable and not limited to a specific family of geometries but rather to all configurations wherein there is enhanced coupling between evanescent modes (near-field) and propagating modes (far-field) -see two specific examples reported recently [21,22]. A recent work by Greffet et al.
[23] critically reviewed publications claiming to exceed Planck's law. Their conclusion was that in the far field regime the widely used term "Super-Planckian" is a misnomer and in cases where such phenomena are reported are in fact cases where the absorption cross section is larger than the geometrical cross section. We shall use the term "Super-Planckian" nonetheless in this paper since it is a widely accepted term in such publications.
In recent years, there has been an increased interest in blackbody radiation mainly due to the emergence of thermal photovoltaic systems (TPV) which are described in great detail in [24][25][26][27][28], where the ability to enhance and control the thermal spectrum translates into higher efficiency of solar cells. TPV systems may benefit immensely from the use of meta-materials [29-31] or structural design of the emitter's surface [32, 33] therefore the ability to investigate with one tool, a large variety of configurations involving blackbody radiation and FDT is evident.
The goal of this study is to present a numerical method designed to calculate thermal radiation produced by a body of arbitrary shape and electromagnetic properties. Several numerical methods for this purpose are described in various publications. These methods include the scattering formulation [ , which is also referred to as the Method of Auxiliary Sources. Section 2 briefly recapitulates the derivation of thermal emission based on FDT and Section 3 demonstrates how the two (blackbody radiation and SMT) can be harnessed together. In Section 4, we assess the numerical error and the method is tested against a configuration that can be solved analytically; simulation results for several other configurations are also presented. In all cases we focus on scenarios when Planck's law is exceeded and/or the spatial resonances of the body are reflected in the thermal radiation.

Essentials of the formulation
The setup of the problem is illustrated in Fig. 1. A homogenous body of arbitrary geometry composed of non-magnetic material with complex permittivity of j ε ε ε ′ ′′ = − is situated in free space. This body assumed to be in thermodynamic equilibrium at temperature T and emits radiation due to random current distributions fluctuating within the body.
Superscripts e and m indicate whether this is an electric or magnetic current density; Conceptually, with this information and the geometry of the body, it is possible to determine the electromagnetic field in terms of a dyadic Green's function which written in operator form reads thus, the Poynting vector associated with the far field can be easily established. Furthermore, since we focus on the far-field component, it is natural to employ the reciprocity theorem [2,12]. Specifically, rather than calculating the field due to "thermal" current densities and with it, the power in the body, we determine the power absorbed by the body when exposed to radiating dipoles located far away from it. Furthermore, we assume an enclosing sphere of radius R ∞ which is much larger than the wavelengths of interest as well as the largest dimensions of the body; its center coincides with that of the body. On this sphere, we place a set of two electric and two magnetic unit dipoles pointing in the θ and φ directions as depicted in Fig. 1.
in order to sample the tangential EM field components using the reciprocity theorem using which the radial Poynting vector and consequently the total radiated power can be found.
Using these dipoles, the tangential components of the thermal EM fields in the body can be determined through the reciprocity theorem in Eq. (3) where the superscript "th" denotes thermal fields or currents originating from the blackbody and the superscript "test" relates to the fields of the test dipoles or their currents.
Using expressions for the cross-correlation of the currents in the body Eq. (1) and Eq. (3), Rytov determined the spectral density of the radiated radial component of the Poynting vector which in our notation has the form Here we employed Eq. (1) and as a result, only the diagonal terms of the dyadic Green's function play a role and therefore we redefined are the "electric" fields generated by the electrical unit-dipole pointing in the θ or φ directions and located on the sphere of radius R λ ∞ >> , at the observation point -in other words these are Green's function components where the source is on the sphere and the observation point is in the body. The next major step is to find a way to numerically determine these Green's operators.

Source model technique (SMT)
Because we are dealing with an arbitrarily-shaped body, the normalized electrical field components in Eq. (4) need to be evaluated numerically. We chose the SMT, which is described in detail in [43]. In the SMT, the fields inside each homogeneous region are approximated by a superposition of fields of fictitious elementary sources whose fields are known analytically (Hertzian dipoles in our case) located outside of this region. Hence, integrations over current distributions typical of other integral methods are avoided and the continuity of the fields across region boundaries can be imposed in the simple point-matching sense. The method is efficient in the sense that it affords a good trade-off between accuracy and computational resources, it is rigorous as it allows increasing the accuracy by gradually increasing computational resources, and it is flexible since it can accommodate general geometries and analysis scenarios. respectively. Together with the incident field, they must satisfy the continuity of the tangential components of the fields across the surface of the body Γ with normal n . We where ( ) v t represents two orthogonal unit vectors tangential to Γ ; the vector ρ to place the sources that generate the external field on a mathematical surface Γ int confined by Γ and the internal field by a series of sources located on a mathematical surface Γ ext that confines Γ . Each group of sources radiates into its relevant region while operating in homogeneous unbounded space characterized by the dielectric coefficient of the region: for the external field the dielectric coefficient is that of the vacuum whereas for the internal field, the dielectric coefficient is that of the body -see Figs. 2(b)-2(c). Finally, we conveniently express the sources in each group in terms of Hertzian dipoles whose field components, when radiating in an unbounded homogenous space, are known analytically. in which  is a generalized voltage vector related to the values of the incident field at the test points and  is a generalized impedance matrix, which is a function of the free-space Green's functions. Once we determine  , we can readily determine the various operators  , which determine the emitted spectrum via Eq. (4).
An important aspect of the SMT is its ability to control the error. It is done by searching for a fictitious source placement configuration which would provide an error below a user predefined threshold. In the framework of this study we used an algorithm which performs such search whilst attempting to minimize the amount of fictitious sources to reduce simulation time.

Results and discussion
To test the validity of SMT for solving blackbody radiation problems we ran a set of simulations for a radiating sphere of radius R and compared it to an analytical solution developed in [44]. The method was found to be suitable for relatively small bodies ( )

.
R λ < For larger dimensions it requires a very large number of fictitious sources to maintain a sufficiently low error. Thus the calculation task becomes more computationally demanding both in terms of memory and runtime (roughly of the order of ( ) 2 O n ). Figure 3 illustrates the error as a function of a sphere's radius to wavelength ratio. It is apparent that even for spheres whose radius to wavelength ratio is ~2 we need a fairly large amount (~1600) of sources if we are to maintain an error below 1%. In all the simulations that follow, the far-field power was computed and maximal error threshold was set to be 0.5%. Figure 4 compares to the analytic results (solid lines) and Planck's law (black line). There is almost a perfect match up between the numerical and theoretical calculations-the total error does not exceed 0.5%. From the physics perspective, it is evident that for large spheres the location of the thermal wavelength (peak of the distribution) and the general shape of the spectrum are becoming similar to Planck's law while for smaller dimensions, we observe greater deviations from it and for certain configurations the Planckian spectrum is exceeded. Moreover, the body behaves as if its temperature is higher as revealed by the curve which represents the 0.1 m μ radius sphere. Figure 4 for several values of the conductivity and they are compared to Planck's law. It is noted that the emitted power strongly depends on the conductivity and for specific values it may exceed Planck's law. It is also interesting to realize that for poorly conducting bodies the spectrum contains peaks which coincide well with the resonant frequencies in a spherical cavity [45] (see also inset). For super-Planckian conditions where the conductivity is high, these resonances are not observed such as in the case of ( ) 2. ε = r Re Fig. 4(c) shows the dependence the spectrum on the real part of the dielectric coefficients. As the dielectric coefficient is elevated, the resonant character of the emission is more pronounced and the excess of power (above Planck's law) is also highly increased. In addition to the dielectric coefficient, we specify the total emitted power normalized to Planck's values: . It is evident that for each radius there is an optimal conductivity for which the total radiated power is maximal. It is also evident for a certain range of radii Planck's law can be exceeded and that for very small bodies and very large bodies it serves at an upper limit.   Next step is to consider the thermal radiation of a conducting sphere with a spherical concentric void. Figures 6(a) and 6(b) show the emitted spectrum and the power -the latter is normalized to a full sphere of the same external radius. The spectrum may exceed both Planck's value as well as that of a full sphere. Figures 6(c) and 6(d) are the same plots for lower conductivity and higher dielectric coefficient. The resonant character of the spectrum is evident in this case too.
Lastly, to demonstrate that thermal radiation, contrary to what is commonly believed, does not depend on global geometrical properties such as surface area or volume but rather on the specific geometry of the body, we analyzed the radiation of ellipsoids of revolution. It is possible to find a pair of ellipsoids with identical volume and surface area but different semiprincipal axes such that one is a prolate and the other is an oblate ellipsoid as shown in Fig.  7(c). Figure 7(a) shows the spectral power density of such two ellipsoids. It is evident that despite them having the same volume and surface area, their emitted power spectra differ greatly one from another and differ from the spectra of spheres with same surface area. It can also be seen that for such miniscule dimensions the calculated results are very different than what Planck's law predicts.
The SMT also allows us to compute the thermal radiation pattern of blackbodies. Figure  7(b) shows the radiation pattern of the oblate and prolate ellipsoids. We can see that the specific geometry of the body affects not only the spectrum but also the directivity of the radiation.
Vol. 25, No. 12 | 12 Jun 2017 | OPTICS EXPRESS A598 To conclude, in this study we have shown that (i) the SMT can be used to numerically solve thermal radiation problems of bodies of dimensions smaller or comparable to the radiated wavelengths. For larger dimensions this method is becoming increasingly slower because a higher number of fictitious sources and test points are required in order to satisfy the boundary conditions. This significantly increases the computer's runtimes and memory consumption. We demonstrated the fact that (ii) the power spectrum of a small body may contain resonant frequencies, however those are observed for poorly conducting spheres where the total emitted power is significantly lower than Planck's law -for low dielectric coefficients. When considering high values of the dielectric coefficient ( )

Re 1
ε >> the resonant character of the thermal radiation is clearly revealed. It has been shown that (iii) blackbody radiation is not solely affected by global geometrical properties such as surface area or volume, but rather by the specific geometry of the body such as in the case of the ellipsoids. Finally, we have also shown that the (iv) total radiated power can be enhanced by creating a hollow cavity inside the radiating body. For a sphere roughly the size of the emitted wavelength, we were able to enhance the radiated power by over 40% by optimizing the dimensions of the cavity.