Numerical methods and approximations for thermal radiation in periodic gratings

From Thermal-FluidsPedia

Revision as of 08:11, 9 March 2012 by Kpark75 (Talk | contribs)
(diff) ← Older revision | Current revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Some commonly used methods for the study of radiative properties of periodic gratings are described in this section. The focus is on the rigorous couple-wave analysis (RCWA) and effective medium theory (EMT).


Rigorous coupled-wave analysis (RCWA)

Originally developed in the 1980’s, RCWA is one of the efficient tools for modeling radiative properties and analyzing diffraction efficiency of periodic gratings bounded by two semi-infinite media [1]. The gratings can be divided into planar thin slabs to approximate the grating profile to an arbitrary degree of accuracy. Due to the periodic structure the dielectric function in the grating region is expanded in a Fourier series. The EM fields are also expressed as a Fourier series expansion in terms of the spatial harmonics, determined by the Bloch-Floquet condition. The accuracy of the solution depends only on the number of terms retained in the Fourier expansion of the dielectric function and the EM fields. By matching boundary conditions for the tangential components of both electric and magnetic fields at the interfaces of each layer, the coupled-wave equations can be formulated and solved utilizing a state-variable method, which has been significantly improved in the more recent papers. RCWA has also been extended to analyze 2D periodic structures.

A 1D binary grating is shown in Fig. 1 to illustrate the RCWA method in solving Maxwell’s equations for an incident plane wave with a wavevector k from Region I. The grating (Region II) is assumed to form on a semi-infinite medium (such as a substrate), Region III. The dielectric functions εA and εB represent the binary grating with a period (Λ). One can treat Region II as an inhomogeneous medium with a periodic dielectric function in the x direction. The grating height (h) is the thickness of Region II. The filling ratio of medium A is given by   f = w / Λ, where w is the width of the grating and usually medium B is either air, vacuum or a dielectric. The direction of incident wavevector can be described by a zenith angle (θ) and an azimuthal angle (ϕ). The plane of incidence is defined by the direction of incidence and the z axis, except at normal incidence when the y-z plane is taken as the plane of incidence. For linearly polarized incident wave, the polarization status is determined by ψ, the angle between the electric field vector and the plane of incidence. The incident electric-field vector E is given by

E=E_{{\rm i}} \exp \left(ik_{x} x+ik_{y} y+ik_{z} z\right)


where Ei is the incident electric field vector at the origin, and the components of the incident wavevector are given by kx=k sinθ cosϕ, ky=k sinθ sinϕ and kz=k cosθ, with k=2π / λ, where λ denotes the vacuum wavelength. In Eq. (1), the time harmonic term exp(iωt), where ω is the angular frequency, is omitted for simplicity. The incident electric field can be normalized so that

E_{{\rm i}} =\left(\cos \psi \cos \theta \cos \phi -\sin \psi \sin \phi \right)\hat{x}  +\left(\cos \psi \cos \theta \sin \phi +\sin \psi \cos \phi \right)\hat{y}-\cos \psi \sin \theta \, \hat{z}


According to the Bloch-Floquet condition [1], the wavevector components of the jth diffraction order in Region I are given by

k_{xj} =\frac{2\pi }{\lambda } \sin \theta \cos \phi +\frac{2\pi }{\Lambda } j


{k_{y} =\frac{2\pi }{\lambda } \sin \theta \sin \phi }


k_{zj}^{{\rm r}} =\left\{\begin{array}{c} {\sqrt{k^{2} -k_{xj}^{2} -k_{y}^{2} } {\rm ,\; }\, \, \, \, \, k^{2} >k_{xj}^{2} +k_{y}^{2} } \\ {i\sqrt{k_{xj}^{2} +k_{y}^{2} -k^{2} } ,\, \, \, \, \, \, k_{xj}^{2} +k_{y}^{2} >k^{2} } \end{array}\right.


As required by the phase-matching condition, the parallel components of wavevector kxj and ky must be the same for the diffracted waves in all three regions. In Eq. (3c), superscript r refers to the diffraction waves reflected back to Region I. Denote the wavevectors for the jth order reflected or transmitted waves as k_{{\rm r}j} =(k_{xj} ,k_{y} ,k_{zj}^{{\rm r}} ) and k_{{\rm t}j} =(k_{xj} ,k_{y} ,k_{zj}^{t} ) respectively. For transmitted diffraction waves, k_{zj}^{r} in Eq. (3c) can be replaced by k_{zj}^{t} after substituting k_{\rm III} =k\sqrt{\varepsilon _{{\rm III}} } for k. From Eq. (3b), the y component of the wavevector is the same for all diffraction orders and in all three Regions. Hence, the wavevectors for all diffracted waves end on the semi-circles, which are intersects of the plane of constant ky and the hemispherical surfaces with a radius kI and kIII in each half plane. Taking z < 0 or Region I for example, the reflected wavevectors form a half-conical surface. Furthermore, the x components of the wavevectors vary by multiples of 2π / Λ according to Eq. (3a). For higher diffraction orders, the z component of the wavevector becomes purely imaginary, such that the diffracted waves are evanescent. If ϕ = 0º or 180º, all the reflected and transmitted diffraction rays lie in the same plane as the plane of incidence. The electric field in region I can be expressed as the sum of the incidence and all reflected waves (including evanescent waves) and that in region III can be expressed as a sum of all transmitted waves (including evanescent waves).

The magnetic fields in Regions I and III can be obtained from Maxwell’s equation,  H\, =(i\omega \mu _{0} )^{-1} \nabla \times E where μ0 is the magnetic permeability of vacuum. In Region II, the electric and magnetic fields can be expressed as a Fourier series. Due to the periodicity, the dielectric function in Region II can also be expressed in a Fourier expansion

\varepsilon _{}^{{\rm ord}} (x)=\varepsilon (x)=\sum _{p}\varepsilon _{p}^{{\rm ord}} \exp \left(i\frac{2p\pi }{\Lambda } x\right)


It is essential to express the inverse of the dielectric function in Region II as a separate Fourier expansion, i.e.,

\varepsilon _{}^{{\rm inv}} (x)=\frac{1}{\varepsilon (x)} =\sum _{p}\varepsilon _{p}^{{\rm inv}} \exp \left(i\frac{2p\pi }{\Lambda } x\right)


where {\varepsilon _{p}^{\rm ord}} and \varepsilon _{p}^{{\rm inv}} are the pth Fourier coefficient for the ordinary and inverse of ε(x) as defined in Eqs. (4a) and (4b), respectively; and in general, \varepsilon _{p}^{\rm inv} \ne {1 / \varepsilon _{p}^{\rm ord}} . In the numerical calculation, an upper limit of p can be set such that p=0,\; \pm 1,\; \pm 2,...\; \pm 2q, which means that there are a total of 4q+1 terms in the Fourier series. Sufficient diffraction orders must be employed in the RCWA calculation to achieve high accuracy [1].

Figure 1. Schematic of a one-dimensional (1D) binary grating with a period Λ, width w, and height h for arbitrary plane wave incidence.

Effective Medium Theory (EMT)

In the mid-infrared (IR) or longer wavelengths, the wavelengths can be much longer than the grating period. In such case, the enhancement and suppression of transmission and absorption can be qualitatively (sometimes quantitatively) explained by the effective medium theory (EMT). In EMT, a complicated structure can be homogenized into an effective medium with a specific material property to reduce the computational time with reasonably approximate results [2,3]. The zeroth-order effective medium expression has been used for the design of periodic gratings with antireflection effects or wavelength-selective radiative properties [1]. In the zeroth-order EMT when the plane of incidence is perpendicular to the gratings, the effective dielectric functions in region II are given as

\varepsilon _{{\rm TE}} =f\varepsilon _{{\rm A}} +(1-f)\varepsilon _{{\rm B}}



\varepsilon _{{\rm TM}} =\left(\frac{f}{\varepsilon _{{\rm A}} } +\frac{1-f}{\varepsilon _{{\rm B}} } \right)^{-1}


for TE and TM waves, respectively. It is important to note that the effective dielectric functions depend on the polarization of the incident wave, so that the effective medium behaves differently for different polarizations. EMT can explain the large transmission found in a wavelength-insensitive manner for TM waves, as well as the nearly-zero transmittance for TE waves in the mid-IR. Figure 2 shows the measured transmittance for different polarizations of a subwavelength Au grating microfabricated on Si substrate. The SEM images are shown as inset. The geometric parameters of the grating are Λ = 1000 nm, w = 600 nm, h = 100 nm, and the thickness of the double-size polished Si substrate is 400 μm. A Fourier-transform infrared (FTIR) spectrometer was used to measure the transmittance for both TE and TM waves at near normal incidence in a setup with two wire-grid polarizers. It can be seen that the RCWA prediction is in good agreement with the spectrometric data. The TE wave spectra monotonically decrease towards the long wavelength, and the transmittance value is less than 0.01 at λ > 3.3 μm. For the TM wave, the transmittance is close to that of plain Si for λ > 8 μm, although the measured values are slightly lower than those predicted by RCWA for 8μm<λ< 11μm. The large transmittance for TM waves toward long wavelengths is consistent with the EMT, but the EMT can significantly underpredict the transmittance for TE waves. At short wavelengths, there is a deep valley in the transmittance spectra near 3.4 μm due to Wood’s anomaly. Another dip in transmittance occurs at about 1.7 μm and is also associated with Wood’s anomaly. For TE waves, Wood’s anomaly is manifested as a discontinuity in the transmittance but the effect is very weak.

Figure 2. Polarization-dependent transmittance of a subwavelength grating: (upper) SEM images of the Au grating on a Si substrate; Measured transmittance of a grating sample and comparison with RCWA model considering all diffraction orders or only the zeroth order for TE waves (Lower a) and for TM waves (Lower b).

Finite-Difference Time-Domain (FDTD)

Other methods, such as the finite-difference time domain (FDTD), have also been used to compute the radiative properties of micro/nanostructures [4]. FDTD is a central difference scheme in both time and space domains with the second-order accuracy and has long been used for radar cross-section calculations. An advantage is that some numerical schemes developed to treat wave phenomena in computational fluid dynamics, such as alternating-direction implicit, can be integrated with FDTD to solve Maxwell’s equations in micro/nanostructures. In recent years, the method is being applied to nano/microscale electro-optical devices at optical frequencies due to its flexibility, robustness, and accuracy. There exist quite a few commercial codes based on the FDTD method.


[1] Lee, B. J., Chen, Y.-B., and Zhang, Z. M., 2008, “Transmission Enhancement through Nanoscale Metallic Slit Arrays from the Visible to Mid-Infrared,” Journal of Computational and Theoretical Nanoscience, 5, pp. 201-213.

[2] Zhang, Z. M., 2007, Nano/Microscale Heat Transfer, McGraw-Hill, New York.

[3] Chen, Y.-B., Zhang, Z. M., and Timans, P. J., 2007, “Radiative Properties of Pattered Wafers with Nanoscale Linewidth,” Journal of Heat Transfer, 129, pp. 79-90.

[4] Fu, K., Chen, Y.-B., Hsu, P.-f., Zhang, Z. M., and Timans, P. J., 2008, “Device Scaling Effect on the Spectral-Directional Absorptance of Wafer’s Front Side,” International Journal of Heat and Mass Transfer, 51, pp. 4911-4925.

[75] Chen, Y.-B., Lee, B. J., and Zhang, Z. M., 2008, “Infrared Radiative Properties of Submicron Metallic Slit Arrays,” Journal of Heat Transfer, 130, 082404