ARTICLE IN PRESS
Journal of Quantitative Spectroscopy & Radiative Transfer 93 (2005) 175–183 www.elsevier.com/locate/jqsrt
Thermal radiation in 1D photonic crystals Arvind Narayanaswamy, Gang Chen1 Department of Mechanical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA. Received 1 February 2004; accepted 1 July 2004
Abstract Optical properties of periodic structures, or photonic crystals, have been studied extensively over the last two decades. While thermal radiation from surface grating structures has also been investigated, radiation inside a photonic crystal has not been investigated completely. We have developed a method to analyze thermal radiation in 1D periodic layered media using a modified dyadic Green’s function method and the fluctuation–dissipation theorem. Using the method, thermal radiation between layers in periodic structures made of ultra-thin metallic films is analyzed. The unusual features of thermal radiation in such structures are explained. r 2004 Published by Elsevier Ltd. Keywords: Thermal radiation; Electromagnetic theory; Dyadic green’s functions; Phased array method; Fluctuationdissipation theorem
Optical properties of photonic crystals (PC) have been a topic of great interest to researchers over the last two decades [1,2]. Even before the seminal papers on PCs, thermal radiation from surface microstructures had been investigated [3–5]. Modification of thermal emission by finite photonic band gap materials has been studied in the past [6]. More recently, PCs have been used to tailor the thermal emission for thermophotovoltaic energy conversion purposes [7–9]. Thermal Corresponding author. Fax: +1-617-253-3484. 1
E-mail addresses:
[email protected] (A. Narayanaswamy),
[email protected] (G. Chen). Also to be corresponded to.
0022-4073/$ - see front matter r 2004 Published by Elsevier Ltd. doi:10.1016/j.jqsrt.2004.08.020
ARTICLE IN PRESS 176 A. Narayanaswamy, G. Chen / Journal of Quantitative Spectroscopy & Radiative Transfer 93 (2005) 175–183
Nomenclature a c ex,ey,ez E Ge Gh 2p_ H i J kx kb kk k K r r0 d3r0 rp TU U0 Ul x^ zl,zl+1 eo e e00 eN g Zl la,lb mo o op Y
arbitrary vector speed of light in vacuum unit vectors along the x, y, and z direction Fourier component of electric field dyadic electric Green’s function dyadic magnetic Green’s function Planck constant Fourier p ffiffiffiffiffiffiffi component of magnetic field 1 source current wave vector component along x-axis Boltzmann constant in-plane wave vector kxex+kyey+kzez,Re(kz)40 kxex+kyey+kzez,Re(kz)40 location at which field is computed source location volume element of the source primitive basis vector transfer matrix of unit cell unit cell where flux is to be computed lth unit cell beyond U0 unit vector perpendicular to wave vector z coordinates of interfaces of lth layer permittivity of free space relative dielectric function imaginary part of e high-frequency dielectric constant frictional damping frequency kzl(zl+1zl) eigenvalues of TU permeability of free space angular frequency (rad s1) plasma frequency _o=ðe_o=kb T 1Þ; Bose–Einstein function
Superscripts i H
pffiffiffiffiffiffiffi 1 homogeneous
ARTICLE IN PRESS A. Narayanaswamy, G. Chen / Journal of Quantitative Spectroscopy & Radiative Transfer 93 (2005) 175–183 177
T d
transpose phase shift angle (0 to 2p) complex conjugate
Subscripts i,j,m,n p s l
1,2,3 polarization (TE or TM) source layer lth layer
radiation inside PC, with the PC acting as emitter, has not received as much attention though. In order to analyze such structures with periods of the order of the characteristic thermal wavelength, we cannot use the classical theory of radiative transfer [10]. Radiative transfer problems in the near field, or taking into account interference effects, can be solved by determining the dyadic Green’s function (GF) [11–13] of Maxwell’s equations and characterizing the thermal source using fluctuation–dissipation (FD) theorem [14–18]. The dyadic GF of periodic structures have been calculated earlier in order to analyze the behavior of properties like local density of states [19,20], which affect the emission of sources embedded in the structure, but never to analyze radiative flux or energy transfer. First, the dyadic GF method and FD method are described. Following that, an analysis of thermal radiation in 1D PC made of ultra-thin metallic films is presented. The 1D periodic structure is shown in Fig. 1. The unit cell consists of an emitter and a layer of vacuum. The emitter is at temperature Th. The method we describe is not restricted to the structure shown in Fig. 1, but can be applied to any 1D PC. The Fourier components of the electrical and magnetic fields can be expressed in terms of the dyadic GF as follows (all materials are assumed to be non-magnetic): R Eðr; oÞ ¼ iomo Ge ðr; r0 ; oÞdJðr0 ; oÞd3 r0 R Hðr; oÞ ¼ Gh ðr; r0 ; oÞdJðr0 ; oÞd3 r0
(1)
where Ge ðr; r0 ; oÞ and Gh ðr; r0 ; oÞ; the electric and magnetic dyadic GF due to a delta source at r0 are related by Gh ðr; r0 ; oÞ ¼ rr Ge ðr; r0 ; oÞ; Jðr0 ; oÞ is the Fourier component of the current due R þ1to thermal fluctuations. Here, the Fourier transform of f ðtÞ is defined as F ðoÞ ¼ 1=2p 1 f ðtÞeiot dt: To compute the radiative flux, we need to determine the Poynting vector at any location. In Cartesian coordinates, the Poynting vector depends on axes and iaj: Since products of the type E i H j ; E j H i ; where i and j denote the D Cartesian E the source in this case is stochastic, we are interested in E i H j ; the brackets denoting a statistical ensemble averaging over all possible configurations of the system. Specifically, D E the Poynting vector in the z direction depends on the term E x H y E y H x : Using Eq. (1), we
ARTICLE IN PRESS 178 A. Narayanaswamy, G. Chen / Journal of Quantitative Spectroscopy & Radiative Transfer 93 (2005) 175–183
Fig. 1. 1D periodic system being investigated. Each unit cell consists of 2 layers—an emitter and vacuum in between. rp is the primitive lattice vector of the 1D periodic structure. The unit cells extend to 1 and þ1 in the z direction. Each planar layer extends to 1 and 1 in the x and y directions, which are in the plane of the layers. For purposes of calculating flux, the emitters in Ul (l p 0) are assumed to be at a constant temperature Th. Those emitters in Ul (l40) are assumed to be at 0 K. Thus, only emitters with hatching contribute to emission, whereas the whole structure participates in producing the photonic crystal effect.
D
E i H j
E
express as D E R R E i ðr1 ; oÞH j ðr1 ; oÞ ¼ iomo V d3 r V d3 r0 f 0 0 G eil ðr1 ; r; oÞGT hmj ðr1 ; r ; oÞ J l ðr; oÞJ m ðr ; oÞ g;
(2)
where repeated indices refer to summation over all possible combinations. To characterize the spectral power density of the thermal source, we make use of the FD theorem [17,18], which expresses the correlations of the microscopic fluctuations of current in terms of a macroscopic parameter, the imaginary part of the dielectric function, as o 00 ðoÞoYðo; T Þ dmn dðr r0 Þ; J m ðr; oÞJ n ðr0 ; oÞ ¼ p
(3)
where 00 ðoÞ is the imaginary part of the dielectric function of the source. With expressions for the dyadic GF and Eq. (3), the expression in Eq. (2) can be computed numerically. The problem of determining the radiative flux is converted to a problem of determining the GF in the given case. For the case of layered media, GF technique has been used extensively [11–13]. The GF for layers which do not contain the delta source, which we shall call the homogeneous part, is
ARTICLE IN PRESS A. Narayanaswamy, G. Chen / Journal of Quantitative Spectroscopy & Radiative Transfer 93 (2005) 175–183 179
given by 8 9 0 Al eikl zl eiðkl :rks r Þ x^ ðkzl Þx^ ðkzs Þþ > > > > > > > > ZZ < = þikl zl iðKl :rks :r0 Þ ^ ^ e e ð Þ x ð k Þþ x k B dk dk i l zl zs x y H 0 Gelp ðr; r Þ ¼ 2 ; (4) 0 > C l eikl zl eiðkl :rKs r Þ x^ ðkzl Þx^ ðkzs Þþ > 8p kzs > > > > > > : ; 0 Dl eþikl zl eiðKl :rKs :r Þ x^ ðkzl Þx^ ðkzs Þ ^ where x^ ðkz Þ ¼ e^ ðk z Þ ¼ ky ex kx ey kr for TE waves and x^ ðkz Þ ¼ hðkz Þ ¼ kz =k kx ex þ ky ey kr þ ez kr k for TM waves [11]. For those layers that contain the dsource, the GF for a particular polarization p is modified to ( ZZ 0 eiks :ðrr Þ x^ ðkzs Þx^ ðkzs Þ; z4z0 dk dk i x y H 0 0 ; (5) Gesp ðr; r Þ ¼ Gesp ðr; r Þ þ 2 0 8p kzs eiKs :ðrr Þ x^ ðkzs Þx^ ðkzs Þ; zoz0 0 where superscript H indicates that the GF is the homogeneous part; and GH esp ðr; r Þ is of the form as in Eq. (4). The corresponding magnetic field GF can be determined using the relation between electric and magnetic GF. With the above expressions, the unknown constants Al,Bl,Cl, and Dl, and hence the GF for a stack of layers bounded by semi-infinite media, can be determined using the transfer matrix method (TMM) [21]. For infinite periodic structures, though, it is not as straightforward. The method used to determine the GF for 1D periodic media is described briefly below. Normally, for periodic media, Bloch’s theorem is used to determine the band structure, transmissivity, reflectivity, and other such properties. In this case Bloch’s theorem cannot be used because of the presence of the delta source, which disrupts the periodicity. Instead, the problem is circumvented by making use of the phased array method [22]. In this method, a whole array of dsources is introduced in addition to the original one, as shown in Fig. 1. The dyadic GF of this modified problem, termed the periodic GF, is determined initially. The important thing is that the d-source in the Nth unit cell to the right (left) of the unit cell in consideration is phase shifted by Nd(-Nd) and all of them have the same arbitrary orientation, a: It can be shown that the field due to this array of d-sources is given in terms of the periodic GF by Z GeP ðr; r0 ; o; dÞdadðr r0 Þ d3 r0 ; (6) EðrÞ ¼ iomo U0
where the integration domain is restricted to the primary unit cell, U0; and subscript P indicates that the GF is a periodic GF. U0 refers to that unit cell where the d-source is at zero phase. By appealing to superposition, the periodic GF can be expressed in terms of the actual GF as GeP ðr; r0 ; o; dÞ ¼
1 X
Ge ðr; r0 þ nrP ; oÞeind ;
(7)
n¼1
1 ) Ge ðr; r þ nrP ; oÞ ¼ 2p 0
Z
2p
GeP ðr; r0 ; o; dÞeind dd; 0
(8)
ARTICLE IN PRESS 180 A. Narayanaswamy, G. Chen / Journal of Quantitative Spectroscopy & Radiative Transfer 93 (2005) 175–183
where rP is the primitive lattice vector of the 1D periodic lattice as shown in Fig. 1. From Eq. (7), it can be shown that GeP ðr; r0 þ rP ; o; dÞ ¼ GeP ðr; r0 ; o; dÞeid :
(9)
Using Eqs. (4), (5) and (9), TMM can be used to determine the periodic GF for the 1D periodic lattice. The actual GF can be determined from the periodic GF by using Eq. (8). The advantage this method offers in determining the GF for a 1D periodic structure is that the integration mentioned in Eq. (8) can be performed easily using contour integration with z ¼ eid being the complex variable. Since the structure is infinite and symmetric about any receiver, the contribution due to the sources on the left of a plane in vacuum should equal the contribution due to those on the right. Hence the sources are separated into sources on the left and sources on the right of the location at which the flux in the +z direction is to be determined. The characteristics of radiative flux or energy density in such structures are by no means intuitive. Analyzing the radiation in an infinite PC made of metallic films will illustrate that. The dielectric function of the metal in the IR can be approximated by the Drude model as ðoÞ ¼ 1
o2p : oðo þ igÞ
(10)
For our calculation, values of 1 ; op ; and g are 5.0, 6.01, and 0.054 eV, respectively. The thickness of the metal film is 10 nm and the vacuum gap is varied between 200 and 800 nm. Since the medium in-between the metal films is non-absorbing, disregarding near-field effects, the flux is a constant irrespective of the location. The ratio of this flux to the normal flux from a planar black body, a nominal ‘‘emissivity’’, is plotted in Fig. 2. The most striking aspect of the figure is the transition from a lower emissivity region in the long wavelength regions to a higher emissivity in the short wavelength regions. The transition from lower emissivity to higher emissivity can be controlled by varying the thickness of the vacuum layer. The spectral characteristics of the flux can be understood by analyzing the eigenvalues of the transfer matrix characterizing the unit cell of the PC. The transfer matrix for the lth layer relates the tangential electric and magnetic fields on one interface of the layer to the fields on the other interface. It can be obtained from the equality of the tangential electric and magnetic fields at any interface. The transfer matrix of the lth layer of the multilayer film can be written as [21] 2 3 0 12 P 3 Plþ1 l ðo=cÞl cos Z i sin Z l l kzl 6 7 @ 6 7 A (11) 4 5¼ 4 5 for TM polarization; kzl i o=c sin Zl cos Zl ð Þl Qlþ1 Ql 2 6 4
Plþ1 Qlþ1
3
0
7 @ 5¼
cos Zl kzl i o=c ð Þ
sin Zl
12 P 3 l i kzl sin Zl 6 7 A4 5 for TE polarization; cos Zl Ql ðo=cÞ
(12)
ARTICLE IN PRESS A. Narayanaswamy, G. Chen / Journal of Quantitative Spectroscopy & Radiative Transfer 93 (2005) 175–183 181
Fig. 2. ‘‘Emissivity’’ of 1D PC for metal thickness of 10 nm and vacuum thickness of 200 and 800 nm.
where pffiffiffiffi l ðAl þ Bl Þ; Ql ¼ o=ckzl pffiffiffi ðAl Bl Þ ð Þ l pffiffiffiffi Rl ¼ l ðC l þ Dl Þ; S l ¼ o=ckzl pffiffiffi ðC l Dl Þ ð Þ l Pl ¼
Pl ¼ ðAl þ Bl Þ; Ql ¼ Rl ¼ ðC l þ Dl Þ; S l ¼
kzl
ðo=cÞ
for TM polarization;
ðAl Bl Þ
kzl ðC ðo=cÞ l
Dl Þ
for TE polarization:
The transfer matrix for the entire unit cell, T U ; is obtained by multiplying the transfer matrices of successive layers of the unit cell. The transfer matrix T U relates the tangential fields at any location in a unit cell to the tangential fields at a corresponding point in an adjacent unit cell. The eigenvalues of T U are solutions of the quadratic equation: l2 TrðT U Þl þ detðT U Þ ¼ 0:
(13)
Since T U is a unitary matrix, la lb ¼ detðT U Þ ¼ 1; where la and lb are solutions to Eq. (13). The presence of the metallic layer implies absorption. Hence the absolute value of one of the eigenvalues has to be less than 1. Let la be the eigenvalue with absolute value less than 1. Referring to Fig. 1, it can be shown that the coefficients of the mth layer (m=1,2) in the lth unit cell (l=N,y,+N), Alm ; Blm ; C lm ; and Dlm are related to A0m ; B0m ; C 0m ; and D0m by 1 1 0 0 A0m Alm C BB C B B lm C jl j B B0m C ¼ l (14) C C: B a B @ C lm A @ C 0m A Dlm
D0m
Eq. (14) can be re-written as Alm ¼ la Aðl1Þm for l40 and Alm ¼ l1 a Aðl1Þm ¼ lb Aðl1Þm for lo0. Similar relations hold for B, C, and D. In that sense, the eigenvalues of T U can be interpreted as ‘‘transmission’’ coefficients. The expressions for the coefficients indicate that the effect of a delta
ARTICLE IN PRESS 182 A. Narayanaswamy, G. Chen / Journal of Quantitative Spectroscopy & Radiative Transfer 93 (2005) 175–183
source in the 0th unit cell decreases as we proceed from the 0th unit cell in either direction. jla j 1 implies that the radiation originating in a unit cell hardly penetrates into the next unit cell. Hence the emission is roughly that from a thin metallic film. On the other hand, when jla j 1; radiation from a unit cell can penetrate into many other unit cells. Alternatively, the radiation in any unit cell is a result of emission from multiple unit cells, and thereby the effect of the PC can be noticed. A plot of la at various wavelengths is shown in Fig. 3a. The metal film is 10 nm and the vacuum gap is 200 nm. It can be seen clearly that for wavelengths in the lower emissivity regions (2.48, 1.24, and 1.13 mm), the value of jla j 1; whereas for wavelengths in the higher emissivity regions (1.03 and 0.83 mm), jla j 1: It can be seen from Eq. (13) that the characteristics of la is determined by TrðT U Þ alone. The solution to Eq. (13) is given by 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi1 !2 1 1 TrðT U Þ 1A: l ¼ @ TrðT U Þ (15) 2 2 # # #TrðT U Þ=2# determines the nature of ð ð Þ Þ Re ð Tr ð T Þ Þ; the value of In the case when Im Tr T U U # # # # #TrðT U Þ=2#41 Similarly, la : When #TrðT U Þ=2#o1; then jla j 1; implying a higher ‘‘emissivity’’. # # implies jla j 1 and hence lower emissivity regions. A plot#of #TrðT U#Þ=2# for normal incidence is shown in Fig. 3b. The wavelength at which the curve of #TrðT U Þ=2# drops below 1 indicates a transition between lower and higher ‘‘emissivity’’. It can be seen from Fig. 3b that the transition wavelengths correlate well with the transition regions in Fig. 2. To summarize, a method to analyze radiative flux in a 1D PC, with the PC acting as emitter has been developed. The method is based on dyadic Green’s function method, phased array technique, and FD theorem. Using the above method, radiative flux in 1D PC of ultra-thin metallic films as emitters with vacuum layers in-between has been analyzed. The radiative flux in the structure shows selective emission behavior. The eigenvalues of the transfer matrix of the unit
Fig. 3. (left) Eigenvalue (smaller in absolute value) of TU at various wavelengths for metal thickness of 10 nm and vacuum gap of 200 nm.The progression of the curves from longer wavelengths wavelengths mimics the # to shorter # progression from lower ‘‘emissivity’’ to higher ‘‘emissivity’’. (right) A plot of #TrðT U Þ=2# as a function of free space wavelength.
ARTICLE IN PRESS A. Narayanaswamy, G. Chen / Journal of Quantitative Spectroscopy & Radiative Transfer 93 (2005) 175–183 183
cell determine the characteristics of the thermal radiation at a given frequency. By controling the extent of the vacuum gap, the position of the transition between lower emission to higher emission can be shifted. This effect indicates that highly selective thermal emitters can be fabricated using ultra-thin metallic films. Acknowledgements The work is supported by an ONR grant N00014-03-1-0835. The authors also acknowledge the support of NSF in providing a travel grant for presenting this work at the Fourth International Symposium on Radiative Transfer in Istanbul, Turkey. References [1] Yablonovitch E. Inhibited spontaneous emission in solid-state physics and electronics. Phys Rev Lett 1987;58(20):2059–62. [2] John S. Strong localization of photons in certain disordered dielectric superlattices. Phys Rev Lett 1987;58(23):2486–9. [3] Hesketh PJ, Gebhart B, Zemel JN. Organ pipe radiant modes of periodic micromachined silicon surfaces. Nature 1986;324(6097):549–51. [4] Hesketh PJ, Zemel JN, Gebhart B. Polarized spectral emittance from periodic micromachined surfaces. I. Doped silicon: the normal direction. Phys Rev B 1988;37(18):10795–802. [5] Hesketh PJ, Zemel JN, Gebhart B. Polarized spectral emittance from periodic micromachined surfaces. II. Doped silicon: angular variation. Phys Rev B 1988;37(18):10803–13. [6] Cornelius CM, Dowling JP. Modification of Planck black-body radiation by photonic band-gap structures. Phys Rev A 1999;59(59):4736–46. [7] Fleming JG, Lin SY, El-Kady I, Biswas R, Ho KM. All-metallic three-dimensional photonic crystals with a large infrared band gap. Nature 2002;417:52–5. [8] Lin SY, Moreno J, Fleming JG. Three-dimensional photonic crystal for thermal photovoltaic power generation. Appl Phys Lett 2003;83(2):380–2. [9] Heinzel A, Boerner V, Gombert A, Blasi B, Wittwer V, Luther J. Radiation filters and emitters for the NIR based on periodically structures metal surfaces. J Mod Opt 2000;47(13):2399–419. [10] Planck M. The theory of heat radiation. New York: Dover; 1991. [11] Kong TJA, Ding KH. Scattering of electromagnetic waves. New York: Wiley; 2000. [12] Tai CT. Dyadic green functions in electromagnetic theory. New York: IEEE Press; 1993. [13] Chew WC. Waves and fields in inhomegeneous media. New York: IEEE Press; 1995. [14] Nyquist, Thermal agitation of electric charge in conductors. Phys Rev 1928; 32(1): 110–3. [15] Callen HB, Welton TA. Irreversibility and generalized noise. Phys Rev 1951;83(1):34–40. [16] Weber J. Fluctuation dissipation theorem. Phys Rev 1956;101(6):1620–6. [17] Rytov SM, Kravtsov YA, Tatarski VI. Principles of statistical radiophysics, vol. 3. Berlin: Springer; 1987. [18] Landau, Lifshitz. Electrodynamics of Continuous Media. Oxford: Pergamon Press; 1969. [19] Martin OJF, Girard C, Smith DR, Schultz S. Generalized field propagator for arbitrary finite-size photonic band gap structures. Phys Rev Lett 1999;82(2):315–8. [20] Moroz A. Minima and maxima of the local density of states for one-dimensional periodic systems. Europhys Lett 1999;46(4):419–24. [21] Heavens OS. Optical properties of thin solid films. New York: Dover; 1991. [22] Caloz C, Skrivervik AK, Gardiol FE. An efficient method to determine Green’s functions of a two-dimensional photonic crystal excited by a line source—the phased-array method. IEEE Trans Microwave Theory Tech 2002;50(5):1380–91.