In this letter, we consider the question of designing insulator/metal thermovoltaic structures with periodically corrugated interfaces that give optimal performance based on the metric of useful power density. Using a Monte Carlo approach in a robust, rapid, and high-accuracy numerical simulation strategy, we have identified such interface shapes. We searched among the class of sinusoids and found that a flat-interface configuration could be significantly improved in transverse magnetic polarization. More specifically, we found that (i) the performance improves with increasing corrugation amplitude (ii) up to a maximum, (iii) the shape of the corrugation is largely irrelevant, and (iv) the period of the corrugation should be chosen in connection to the bandgap energy of the photovoltaic cell. For the latter, we provide a simple expression as a starting point for practitioners interested in fabricating such structures.

A thermophotovoltaic (TPV) emitter is a structure that, when heated to an appropriate temperature, emits photons with energies suitable for the generation of electricity by a photovoltaic (PV) cell.1 If solar energy is used to heat the emitter, then the full system is a solar thermophotovoltaic (STPV) system,1,2 although applications of TPVs involving nonsolar-based energy (e.g., waste heat) are also of great interest. Ideally, the emitted photons should have energies near and above the bandgap energy of the relevant PV cell. If the emitter were simply a blackbody emitter, then only the operating temperature would be of relevance, but layered and nanostructured TPV structures can lead to emission properties more tailored for efficient energy generation.3,4

Inspired by Jeon et al.,5 we consider a simple Bragg reflector-tungsten TPV system (see Fig. 1) and use a highly efficient computational electrodynamics procedure to predict optimal surface structuring. By corrugating the interfaces in the emitter, we can enhance the emissivity in the relevant, sub-bandgap range of energies. We follow Ref. 6 where the authors reasoned that such corrugations act like a “graded” material and allow one to modify the emissivity. However, rather than considering their sharply varying sawtooth profiles,6 we focus on smooth periodic profiles to corrugate. Not only is this feature advantageous for our numerical algorithm, but also the response is similar to sawtooth shapes readily generated by modern fabrication techniques.

FIG. 1.

Depiction of the full STPV system (left) with special emphasis on the TPV emitter structure (right).

FIG. 1.

Depiction of the full STPV system (left) with special emphasis on the TPV emitter structure (right).

Close modal

The numerical approach we utilize is a high-order perturbation of surfaces (HOPS) algorithm,7–10 which utilizes surface unknowns, giving it an order-of-magnitude advantage in terms of storage and execution time over volumetric approaches such as finite difference,11 finite element,12 spectral element,13 and spectral14 methods. Furthermore, because of its perturbative character and expression in terms of periodic eigenfunctions of the Laplacian, it has advantages over integral equation approaches:15 there is no need for specialized quadratures, periodization strategies, or iteration schemes for solving dense, nonsymmetric positive-definite systems of linear equations.16,17

From our simulations, we have made a number of important discoveries. First, the introduction of periodic corrugations can significantly increase the useful power density (PD) (4.2651 W/cm2 versus 3.4443 W/cm2) of solar cells in transverse magnetic (TM) polarization. Second, the enhancement grows with increasing corrugation amplitude up to a maximal value beyond which the beneficial effects dissipate. Third, the shape of the corrugation does not appear to be crucial since a simple sinusoid produces roughly the same results as those produced by a sawtooth profile. Finally, our results improve as the period, d, of the corrugation is adapted to the bandgap wavelength, λBG, of the underlying PV cell. Based upon the careful study of the emissivity of our structures, we have derived the rule d λBG as a useful starting point for practitioners building such devices. We note that our results depend strongly on polarization: The difference between flat and corrugated interfaces is minimal in transverse electric (TE) polarization; however, it is sizable in TM.

Figure 2 displays the geometry of the configuration we consider: a y-invariant (M +1)-layered insulator-metal structure. An insulator (vacuum with refractive index nvac = 1) occupies the domain above the uppermost interface {z > g1(x)} and the metal (tungsten) fills {z < gM(x)}. A tungsten-alumina alloy spacer (Al2O3 with nal = 1.7682 with tungsten volume fraction 0 ≤ fW ≤ 1) occupies the second layer, {g2(x) < z < g1(x)}, while a Bragg reflector composed of alternating layers of SiO2 (nsi = 1.4585) and TiO2 (nti = 2.6142) occupies the middle of the structure {gM(x) < z < g2(x)}. We focus upon d-periodic grating interfaces, typically gj(x)=g¯j+asin(2πx/d), and work with triply layered Bragg structures so that M = 5. The structure is illuminated from above by monochromatic plane-wave incident radiation of frequency ω, aligned with the grooves. Factoring out the common term exp(iωt), we choose the reduced electric and magnetic fields as unknowns, which, like the incident field, are quasiperiodic.

FIG. 2.

Bragg-tungsten structure with sinusoidal periodic interfaces.

FIG. 2.

Bragg-tungsten structure with sinusoidal periodic interfaces.

Close modal

In this two-dimensional setting, the time-harmonic Maxwell equations decouple into (M +1)-many scalar Helmholtz problems that govern the TE and TM polarizations.18 We denote the invariant (y) directions of the scattered (electric or magnetic) fields by vm(x, z) in layer m and the incident radiation in the upper layer by vi. We seek outgoing quasiperiodic solutions of

(1)
(2)
(3)

where δm, is the Kronecker delta, k(m) = n(m)ω/c, Nm=(xgm,1)T,τm2=1 in TE, and τm2=(k(m)/k(m+1))2=(n(m)/n(m+1))2,inTM.

The Rayleigh expansions18,19 state, for z>|g1|L,

where αp:=α+2πp/d, γp(0)=(k(0))2αp2, with Im(γp(0))0, and the “propagating modes” are U(0)={pZ|αp2<(k(0))2}. The emissivity (equal to the absorbance) is defined as

We follow a HOPS methodology to solve (1)(3), which successively corrects the flat-interface (Fresnel) approximation.8,16,20–22 The approach begins with the assumption that the shapes of the interface deformations gm satisfy gm = εfm (ε ≪ 1) with fm being sufficiently smooth. The smallness assumption on ε can be removed by analytic continuation23,24 and numerically implemented via Padé summation.8,25,26

We utilize the Transformed Field Expansion approach8,20 which we recall here. A simple change in variables has an impact in each layer which flattens the domain interfaces, {z = gj(x)} to {z=g¯j}. This delivers equations for the transformed fields, um(x,z)=vm(x(x,z),z(x,z)), e.g.,

for (1), where forms for Fm can be readily derived.9,10 Classical results7,20 imply that the transformed fields depend analytically on the deformation size ε so that, e.g.,

(4)

This expansion is inserted into the transformed version of (1)(3), resulting in a system of coupled inhomogeneous Helmholtz problems governing um,n to be solved. For the purpose of approximating these, we have selected a stable, High-Order Spectral Fourier-Legendre Galerkin approach14 where

and L are appropriately scaled and translated Legendre polynomials. For smooth profiles fm, the scattered fields, um, are jointly analytic in x, z, and ε, so the coefficients ũm,n,p, decay exponentially fast as {m, n, p, } grow.8 Thus, only a small number of these are required to deliver a high-fidelity solution which can be discovered very quickly. Naturally, as the deformation size ε becomes larger, more perturbation orders n are required, which slows our algorithm.

A figure of merit that measures the utility of our structures comes from Planck's law for the spectral radiance of a black body at temperature T

where is Planck's constant, c is the speed of light in a vacuum, and kB is the Boltzmann constant. With this, the (useful) power density is given by

(5)

where we restricted the integration domain to be {λmin < λ < λmax}. Jeon et al.5 introduced a second figure of merit, the “spectral efficiency,” which, in our simulations, decreased only slightly with the introduction of corrugations.

In all of our simulations, we chose the physical parameters λBG = 2.254 μm, λmin = 0.6 μm, λmax = 6.0 μm, T = 1700 K, and fW = 0.75, where we have used a Maxwell formula to estimate the permittivity of the alloy.5 Geometrically, we have chosen the base layer thicknesses to be 20 nm for tungsten-alumina alloy, 255 nm for SiO2, and 150 nm for TiO2.

To begin, we consider a selection of representative configurations with sinusoidal corrugations of the spatial period, d = 2.5 μm, and amplitudes a = 0, 0.15, 0.275, and 0.35 μm. The results are summarized in Table I. The final column, δ, is a dimensionless measure of the energy defect when the tungsten and alloy layers are replaced by a dielectric. If δ = 0, then energy is perfectly conserved (as it should be in a dielectric structure), while the values of 10−2 to 10−3 indicate that 2–3 digits can be trusted. From this table, we learn a number of things. First, we can increase the useful power density by nearly 25% with the introduction of periodic sinusoidal corrugations. To see why this is the case, we display, in Fig. 3, the emissivity, ϵs, for flat (a = 0 μm) and corrugated (a = 0.25 μm) interface configurations. We notice the significant enhancement of the magnitude of the emissivity below the bandgap wavelength in the presence of corrugations. In addition, from Table I, we note that the improvement in useful power density increases as the corrugation amplitude increases to a maximum around a 0.275 μm but then decreases again as a increases.

TABLE I.

Numerical results for six-layer structures.

a (μm)PD (W/cm2)ΔNxNzN
3.4443 4 × 10−16 20 20 20 
0.15 3.8381 3 × 10−3 20 20 20 
0.275 4.2651 8 × 10−3 20 40 40 
0.35 4.1651 1 × 10−2 20 40 40 
a (μm)PD (W/cm2)ΔNxNzN
3.4443 4 × 10−16 20 20 20 
0.15 3.8381 3 × 10−3 20 20 20 
0.275 4.2651 8 × 10−3 20 40 40 
0.35 4.1651 1 × 10−2 20 40 40 
FIG. 3.

Plot of the emissivity, ϵs, versus wavelength, λ, for the sinusoidally corrugated interface configuration with a = 0 and 0.25 μm. (Nx = Nz = 20 and N = 20).

FIG. 3.

Plot of the emissivity, ϵs, versus wavelength, λ, for the sinusoidally corrugated interface configuration with a = 0 and 0.25 μm. (Nx = Nz = 20 and N = 20).

Close modal

To discern the effect of the particular shape of the interface deformations, we ran these simulations again with a = 0.18 μm and interfaces shaped as sawtooth profiles (see Fig. 4). These simulations revealed that the useful power density increases from 4.01 W/cm2 to 4.03 W/cm2 when the sharp sawtooth corrugations are introduced, which is clearly not significant.

FIG. 4.

Bragg-tungsten structure with sawtooth periodic interfaces.

FIG. 4.

Bragg-tungsten structure with sawtooth periodic interfaces.

Close modal

Continuing, we decided to investigate the possibility of improving the performance of our design by varying the period of our corrugations. For this, we selected 100 samples of d from the uniform distribution U(1.6, 3.0) (in microns) for a fixed amplitude a = 0.18 μm; the results are depicted in Fig. 5. Here, we see a significant spike near the value d = 2.8 μm. Our explanation for this result is that it is within a small neighborhood of λBG so that our emissivity profile is tailored to the bandgap of the PV. To make this more quantitative, we display the emissivity for a period d near this critical value λBG in Fig. 6 and approximately below/above this value in Fig. 7 with amplitude a = 0.18 μm corrugations. Here, we see how the emissivity curve is “optimal” in the d = 2.88 ≈ λBG case with all of the plasmonic peaks located at wavelengths up to, but not past, the bandgap wavelength. By contrast, for d <λBG, the peaks do not exist all the way up to the bandgap wavelength, while for d >CλBG, there are “wasted” peaks beyond.

FIG. 5.

Useful power density versus period of the sinusoidally corrugated interfaces (Nx = Nz = 20 and N = 20).

FIG. 5.

Useful power density versus period of the sinusoidally corrugated interfaces (Nx = Nz = 20 and N = 20).

Close modal
FIG. 6.

Emissivity versus wavelength for period d = 2.88 μm (d λBG) corrugated interfaces (Nx = Nz = 20 and N = 20).

FIG. 6.

Emissivity versus wavelength for period d = 2.88 μm (d λBG) corrugated interfaces (Nx = Nz = 20 and N = 20).

Close modal
FIG. 7.

Emissivity versus wavelength for period d = 2.4 μm (d <CλBG) and d = 3.2 μm (d >CλBG) corrugated interfaces (Nx = Nz = 20 and N = 20).

FIG. 7.

Emissivity versus wavelength for period d = 2.4 μm (d <CλBG) and d = 3.2 μm (d >CλBG) corrugated interfaces (Nx = Nz = 20 and N = 20).

Close modal

In this letter, we addressed the problem of designing TPV emitters consisting of Bragg reflectors overlaying tungsten with periodically corrugated interfaces. Using a rapid and robust HOPS solver, we have found such structures using the useful power density figure of merit. We determined that corrugations of simple sinusoidal form can give significant enhancements which grow with increasing amplitude up to a maximal value. The shape of the corrugation did not appear to be particularly important; however, the period, d, should be adapted to the bandgap wavelength of the PV cell, λBG, with a useful first approximation expressed by d λBG.

We note that our calculations focus on normal incidence (emission), which is common in the literature. As discussed in the supplementary material of Ref. 5, although angles near normal are expected to be most important, a more complete calculation would involve integrating over all the angles of emission. The conditions for plasmonic/diffractive resonances depend on the angle [e.g., Eqs. (1) and (2) of Ref. 27] so that there will likely be some deterioration in the overall power densities.

This paper has focused on achieving emissivities that can enhance power density associated with a PV cell at a particular band-gap energy. However, another important aspect of solar energy related problems (e.g., concentrating solar power28) is to design structures that can absorb over the entire solar spectrum. Figure 7 shows that one can achieve absorption (emissivity) at somewhat longer wavelengths by increasing the periodicity. The extent to which one could optimize the periodicity represents an interesting problem that we plan to address.

We thank J. Foley for helpful discussions. This material is based upon the work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, under Contract No. DE-AC02-06CH11357, and the National Science Foundation through Grant No. DMS-1522548 (D.P.N.). The use of the Center for Nanoscale Materials, an Office of Science user facility, was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357.

1.
P.
Bermel
,
M.
Ghebrebrhan
,
W.
Chan
,
Y. X.
Yeng
,
M.
Araghchini
,
R.
Hamam
,
C. H.
Marton
,
K. F.
Jensen
,
M.
Soljacic
,
J. D.
Joannopoulos
,
S. G.
Johnson
, and
I.
Celanovic
,
Opt. Express
18
,
A314
(
2010
).
3.
C.
Ungaro
,
S. K.
Gray
, and
M. C.
Gupta
,
Appl. Phys. Lett.
103
,
071105
(
2013
).
4.
J. J.
Foley
,
C.
Ungaro
,
M. C.
Gupta
, and
S. K.
Gray
,
Opt. Express
23
,
A1373
(
2015
).
5.
N.
Jeon
,
J.
Hernandez
,
D.
Rosenmann
,
S.
Gray
,
A.
Martinson
, and
J.
Foley
,
Adv. Energy Mater.
8
,
1801035
(
2018
).
6.
C.
Ungaro
,
S.
Gray
, and
M.
Gupta
,
Opt. Lett.
39
,
5259
(
2014
).
7.
D. P.
Nicholls
and
F.
Reitich
,
Proc. R. Soc. Edinburgh, Sect. A
131
,
1411
(
2001
).
8.
D. P.
Nicholls
and
F.
Reitich
,
J. Opt. Soc. Am. A
21
,
606
(
2004
).
9.
Y.
Hong
and
D. P.
Nicholls
,
J. Comput. Phys.
330
,
1043
(
2017
).
10.
Y.
Hong
and
D. P.
Nicholls
,
J. Comput. Phys.
345
,
162
(
2017
).
11.
R. J.
LeVeque
,
Finite Difference Methods for Ordinary and Partial Differential Equations
: Steady-State and Time-Dependent Problems (
Society for Industrial and Applied Mathematics (SIAM)
,
Philadelphia, PA
,
2007
), pp.
xvi+341
.
12.
C.
Johnson
,
Numerical Solution of Partial Differential Equations by the Finite Element Method
(
Cambridge University Press
,
Cambridge
,
1987
), p.
278
.
13.
M. O.
Deville
,
P. F.
Fischer
, and
E. H.
Mund
,
High-Order Methods for Incompressible Fluid Flow
, Cambridge Monographs on Applied and Computational Mathematics Vol.
9
(
Cambridge University Press
,
Cambridge
,
2002
), pp.
xxviii+499
.
14.
J.
Shen
,
T.
Tang
, and
L.-L.
Wang
,
Spectral Methods
, Springer Series in Computational Mathematics, Algorithms, Analysis and Applications Vol.
41
(
Springer
,
Heidelberg
,
2011
), pp.
xvi+470
.
15.
D.
Colton
and
R.
Kress
,
Inverse Acoustic and Electromagnetic Scattering Theory
, 2nd ed. (
Springer-Verlag
,
Berlin
,
1998
), pp.
xii+334
.
16.
D. P.
Nicholls
,
J. Opt. Soc. Am., A
32
,
701
(
2015
).
17.
D. P.
Nicholls
,
SIAM J. Numer. Anal.
55
,
144
(
2017
).
18.
R.
Petit
, ed.,
Electromagnetic Theory of Gratings
(
Springer-Verlag
,
Berlin
,
1980
), pp.
xv+284
.
19.
P.
Yeh
,
Optical Waves in Layered Media
(
Wiley-Interscience
,
2005
), Vol.
61
.
20.
D. P.
Nicholls
and
F.
Reitich
,
J. Opt. Soc. Am. A
21
,
590
(
2004
).
21.
D. P.
Nicholls
,
F.
Reitich
,
T. W.
Johnson
, and
S.-H.
Oh
,
J. Opt. Soc. Am., A
31
,
1820
(
2014
).
22.
D. P.
Nicholls
,
S.-H.
Oh
,
T. W.
Johnson
, and
F.
Reitich
,
J. Opt. Soc. Am., A
33
,
276
(
2016
).
23.
D. P.
Nicholls
and
F.
Reitich
,
Numer. Math.
94
,
107
(
2003
).
24.
B.
Hu
and
D. P.
Nicholls
,
Proc. R. Soc. Edinburgh A
140
,
367
(
2010
).
25.
O.
Bruno
and
F.
Reitich
,
J. Opt. Soc. Am. A
10
,
2307
(
1993
).
26.
D. P.
Nicholls
and
F.
Reitich
,
J. Comput. Phys.
170
,
276
(
2001
).
27.
H.
Gao
,
J. M.
McMahon
,
M. H.
Lee
,
J.
Henzie
,
S. K.
Gray
,
G. C.
Schatz
, and
T. W.
Odom
,
Opt. Express
17
,
2334
(
2009
).