Near-field thermophotovoltaics for efficient heat to electricity conversion at high power density

Thermophotovoltaic approaches that take advantage of near-field evanescent modes are being actively explored due to their potential for high-power density and high-efficiency energy conversion. However, progress towards functional near-field thermophotovoltaic devices has been limited by challenges in creating thermally robust planar emitters and photovoltaic cells designed for near-field thermal radiation. Here, we demonstrate record power densities of ~5 kW/m2 at an efficiency of 6.8%, where the efficiency of the system is defined as the ratio of the electrical power output of the PV cell to the radiative heat transfer from the emitter to the PV cell. This was accomplished by developing novel emitter devices that can sustain temperatures as high as 1270 K and positioning them into the near-field (<100 nm) of custom-fabricated InGaAs-based thin film photovoltaic cells. In addition to demonstrating efficient heat-to-electricity conversion at high power density, we report the performance of thermophotovoltaic devices across a range of emitter temperatures (~800 K–1270 K) and gap sizes (70 nm–7 µm). The methods and insights achieved in this work represent a critical step towards understanding the fundamental principles of harvesting thermal energy in the near-field.


Supplementary Note 1. Fabrication of the suspended thermal emitter
A schematic diagram of the fabrication process for the emitter device is shown in Supplementary   Fig. 1a. (Step 1) A commercially available 4" SOI wafer (Ultrasil LLC) with a heavily doped (Ptype B-doped) device layer and a thickness of 60 ± 1 μm is chosen. The specified resistivity of the top device layer is <0.005 Ω.cm, which corresponds to a doping level ~3 × 10 19 cm -3 . The handle layer is undoped with a thickness of 500 ± 10 μm. The oxide layer thickness between the device layer and the handle layer is 2 μm ± 5%, which acts as an etch stop for fabricating the suspended emitter. (Step 2) The top surface is patterned using standard lithographic techniques to define the mesa surface. This pattern is etched to a depth of 15 μm in a deep reactive ion etching (DRIE) tool. Next, we proceed to release the emitter device to form a suspended structure. The frontside of the wafer is first coated with a thick photoresist (AZ-9260) for protection. The backside is subsequently coated with a thick photoresist (10 μm) and patterned to expose windows for etching.
The exposed silicon undergoes a through-etch (DRIE -500 μm deep) process until the stop BOX layer is reached. The residual BOX layer is then etched using a recipe that selectively etches silicon dioxide, resulting in a suspended monolithic silicon emitter. Next, the released chips are thoroughly cleaned in hot Piranha to remove any organic contaminants while retaining the smoothness of the mesa surface. (Step 5) A 10 nm-thick layer of aluminum nitride (AlN) is coated using an atomic layer deposition (ALD) tool. In Supplementary Fig. 1b we show a SEM image of a fabricated silicon emitter device, suspended from the base through two silicon beams.

Supplementary Note 2. Fabrication of the InGaAs photovoltaic cell
A schematic diagram of the fabrication process for the PV cell is shown in Supplementary Fig. 2.

Thermal characteristics
To understand the thermal and mechanical characteristics of the emitter at high temperatures, we performed FEM analysis, using COMSOL Multiphysics, on the structure shown in Supplementary   Fig. 3a. Specifically, we employed the heat transfer, solid mechanics, electric currents and surfaceto-surface radiation modules, which were coupled to the physics modules of thermal expansion and electromagnetic heating by using the relevant boundary conditions. In the heat transfer module, a constant temperature (300 K) condition was applied on surfaces A and B (see Supplementary Fig. 3a). Bipolar voltages of V+(~3.75) and V-(~-3.75 V) were also applied on surfaces A and B respectively, which were then used to solve for the resulting currents to determine the Joule heating in the electromagnetic heating module. To account for the radiative coupling to the environment, a nominal emissivity of 0.7 is applied on all the surfaces of the emitter, while the environmental temperature was set to 300 K. The thermal conductivity of silicon as a function of temperature was approximated from that of bulk intrinsic silicon as described in ref. 1 . Supplementary Fig. 3b shows the resulting temperature rise (color scale in kelvin) in the emitter structure, with the mesa being at the highest temperature of 1292 K. The distributed Joule heating in the beams results in a temperature gradient that is nearly exclusively spread along the length of the beams, while the mesa surface is seen to be isothermal with less than 4 K temperature change along the two lines shown in Supplementary Fig. 3b. A transient thermal analysis was performed on the structure to estimate the time constant of the emitter. The temperature change with respect to time is shown in Supplementary Fig. 3c, indicating that the thermal time constant is ~20 ms.
We note that all the near-field measurements were carried out after allowing for a thermal stabilization time of at least 3 s.
An important consideration in the design of thermal emitters is the possible distortion of the emitter due to thermal expansion and bi-material effects. We note that by suspending a monolithic silicon emitter, instead of creating a doubly clamped structure; we reduce potential bending and buckling effects. In Supplementary Fig. 3d, we show the computed absolute displacements for an emitter heated to 1292 K. Thermal expansion occurs predominantly in the direction along the length of the beams with a maximum displacement of ~1 μm at the distal end of the cantilever (Supplementary Fig. 3d) and the mesa surface is expected to remain flat at all temperatures. We note that, for simplifying the simulations, we did not consider the temperaturedependence of resistivity and emissivity of doped silicon in this simulation, but this simplification is not expected to affect any of our conclusions.
Finally, we note that the heat dissipated in the beams is lost through conduction along the beams and radiation from all the surfaces to the surroundings. The contribution of different pathways is affected by the temperature and the gap size. For example, at Temitter = 1270 K achieved by dissipating 550 mW in the emitter and a gap size of 100 nm, the total radiative heat from the surfaces facing the PV cell is 2 mW (0.4%) and from the backside and the side walls of the emitter is ~10 mW (1.8%). The major heat pathway is through conductance along the beams of 537 mW (98%).

Mechanical stiffness
To precisely maintain the gap size between large planar surfaces, it is important to design a stiff cantilever to minimize deflections at the smallest gap sizes and avoid snap-in (for example due to electrostatic attraction) at nanometer-scale gap sizes. Our emitter consists of a large circular island suspended by two beams, each 20 µm wide, 270 µm long and 45 µm deep. The resulting stiffness of the emitter as a function of temperature was also calculated using a FEM model similar to that described in the previous section. The temperature-dependent elastic modulus (E) of silicon is estimated from the semi-empirical formula obtained from ref. 2 Here, T is the temperature of silicon, E0 = 167.5 GPa is the elastic modulus at 0 K, B = 15.8 MPa K -1 and T0 = 317 K are the fitting parameters 2 . A constant force, F = 10 mN, was applied on mesa normal to the surface ( Supplementary Fig. 3a), while the emitter is heated, using the model described in the previous section, by applying a voltage V that is varied from 1.5 V to 3.9 V in Fig. 3e) is monitored at these temperatures and the stiffness is calculated as k = F/ δ ( Supplementary Fig. 3f). This was found to be ~2.15 kN/m at 800 K and reduced by ~5% as the temperature is increased to 1300 K.

Supplementary Note 4. Electro-thermal characteristics of emitter
In this study, controlled heating of the emitter is achieved by employing heavily-doped silicon (~3 × 10 19 cm -3 ) as the resistive material. The electrical circuit used in this study involves applying a bipolar voltage using a DC power supply and measuring the electrical current through the circuit using a multimeter, as shown in the inset of Supplementary Fig. 4a. In Supplementary Fig. 4a, we plot the resistance of the emitter as a function of the power dissipated, while the corresponding temperature of the emitter is plotted on the secondary axis (the details of this temperature measurement are provided in supplementary section 6). An initial increase in the resistance of the doped silicon is observed when the emitter is heated, with the resistance reaching a maximum at ~1200 K, beyond which it starts to decline. This trend can be understood by considering the temperature-dependent electrical resistivity (ρ(T)) of doped silicon as , where n(T) is the temperature-dependent number of charge carriers, e is the elementary charge and μ(T) is the temperature-dependent charge carrier mobility. At relatively low temperatures, ρ(T) is dominated by the temperature dependency of carrier mobility (μ(T)) which decreases with increasing temperature 3 . As the temperature is increased, the number of thermally generated intrinsic carriers (n(T)) increases significantly 4 , causing ρ(T) to drop with temperature. Such a behavior, called thermal runaway, is observed in doped silicon resistors as described in ref. 4 . An additional observation is that the resistance curve displays hysteresis after few hours of operation at high temperatures, and thus, it is not feasible to use the resistance as an indicator of temperature, instead, we measure the temperature in an independent experiment (supplementary note 6).
To reduce the electrostatic forces between the parallel plates of the emitter and PV cell during the near-field approach, we apply a bipolar voltage across the device bringing the mesa potential close to ground. To test the effectiveness of this approach, we fabricated a few 3-terminal devices and tested them in representative conditions. In Supplementary Fig. 4b two examples are shown, one by applying a bipolar voltage (case a) and the other by passing current (case b) from a Keithley 6220 current source (insets of Supplementary Fig. 4a) through the two outer terminals of the device, while monitoring the voltage (VMon) on the mesa with the third terminal in the center.
It can be seen that for case a, VMon remains close to zero potential while a continuous increase is seen for case b, indicating that the higher electrostatic force for case b could lead to snap-in of the devices at higher gap sizes than that for case a.

Supplementary Note 6. Emitter Temperature Measurement via Scanning Thermal Microscopy
The temperature of the emitter as a function of dissipated power was measured on the same device after the near-field experiments were performed, using a quantitative high-temperature scanning thermal microscopy (SThM) technique in an ultra-high vacuum environment following an approach described in ref. 5 . Briefly, a custom-fabricated SThM probe with an integrated Pt serpentine resistor that serves both as a heater and thermometer was used for these measurements.
The temperature at the probe's tip (Ttip) can be accurately quantified with the Pt thermometer.
( ) into the tip. In our measurements the probe was placed close to the center of the emitter's mesa but at a distance of ~5 µm above the surface, with the orientation shown in Supplementary   Fig. 6a. Then the probe was moved towards the mesa to contact the emitter. The thermal behavior of the resulting structure can be represented by the DC equivalent thermal circuit shown in Supplementary Fig. 6b, where RC and RTS are the cantilever thermal resistance and tip-sample thermal contact resistance, respectively. Given the value of RTS, the surface temperature (Ts) can be calculated based on the measurement of Ttip using where Ttip-retract and Ttip-contact refer to SThM tip temperatures when retracted (out of contact) and when in contact with the mesa. We note that upon contact of the probe and emitter, an additional heat transfer path is provided by which the emitter's temperature is reduced. However, this decrease is less than 0.4% of the total temperature rise of the emitter due to the significantly larger thermal resistance of the additional path (RTS + RC) compared to the thermal resistance of the emitter. Hence this effect can be neglected. in the AC temperature amplitude when the probe is in contact compared to out of contact provides direct information about the value of RTS based on a single-state dynamic thermal model as described in ref. 5 . Supplementary Fig. 6c shows the histogram of two measurements for an emitter that dissipated powers of 0.314 W and 0.45 W, showing mean temperatures of 414.8 °C and 794.1 °C respectively. As can be seen, the measurement uncertainty increases at higher temperatures, because it is more difficult to achieve a stable contact at elevated temperatures. Finally, the mean value of emitter surface temperature as a function of emitter power is shown in Supplementary   Fig. 6d. Since the temperature of the emitter as a function of dissipated power depends primarily on the temperature-dependent thermal conductivity of silicon, this relationship is expected to hold in our NF-TPV and the SThM temperature measurements. The non-monotonic increase in temperature with power may be attributed to increasing radiative coupling to the environment at higher temperatures and due to the temperature and position-dependent changes in material properties such as resistivity along the beams.

Supplementary Note 7. Modelling of the NF-TPV system
As described in the methods section of the main manuscript, our theoretical modelling involves calculating the radiative energy transfer between the emitter and the PV cell, and further estimating PMPP from a PV cell model. We begin by describing our approach to modelling the spectral energy transfer.

Modelling Radiative Energy Transfer
The geometry and composition of the emitter and the PV cell used for modelling the radiative energy transfer are shown in Supplementary Fig. 7a. First, we model the permittivity of n-InGaAs by combining the intrinsic optical response of pure InGaAs, and free carrier absorption due to doping as: , when , when , where is the angular frequency corresponding to the bandgap of InGaAs.
For , at wavelengths shorter than 12 µm, we use tabulated values 6 . At wavelengths longer than 12 µm, the intrinsic optical response is due to phonons. We model at these long wavelengths by interpolating between the permittivity of InAs 7 and the permittivity of GaAs 7 using the composition (In0.53Ga0.47As): 0.53 0.47 , when 12 .
We use the following Drude term to account for the free-carrier absorption in the n-doped material: where the electron concentration is given by , where the electron mobility μ is obtained from ref. 10 . For a n-doping level 1 10 , the effective mass is * 0.041 , where is the mass of electron, and the electron mobility is 8727 ⋅ ; At an n-doping of 1 10 , the effective mass * 0.0483 , and electron mobility is 5639 ⋅ .
To model the permittivity of doped InP, we use a similar approach by accounting for both the intrinsic optical response and free carrier effects. The permittivity of doped InP is: , when , when , where is the angular frequency corresponding to the bandgap of InP. The permittivity for intrinsic InP is taken from ref. 7 . We use a Drude term to account for the free carrier optical response Here, the electron and hole concentrations are given by where and are donor and acceptor concentrations, respectively, and the intrinsic carrier concentration for InP 8 at room temperature is 1. 3 10 . The electron effective mass * is determined by using the dependence of electron effective mass on electron concentration 11 .
The electron scattering rate is Γ * , where the electron mobility μ is obtained from ref. 12 .
For the hole effective mass (ref. 9) we use * 0.6 . The hole scattering rate is Γ * , where the hole mobility μ 150 ⋅ / and was obtained from ref. 13 .
Based on the above parameters, we adopt the following values for our model: For the ndoped InP layer, with 1 10 and 0, the effective electron mass is * 0.0831 , the electron mobility is 2096 ⋅ , the effective hole mass is * 0.6 , and the effective hole mobility is 150 ⋅ . For the p-doped InP layer, with 1 10 and 0, the effective electron mass is * 0.08 , electron mobility is 5400 ⋅ , the effective hole mass is * 0.6 , and the effective hole mobility is 46.4 ⋅ . Finally, the permittivity for gold is obtained from ref 7 .
Next, we describe the modelling of the emitter, where the circular mesa consists of a 60 μm thick layer of p-doped Si (doping concentration 3×10 19 cm -3 ), covered by 10 nm-thick layer of AlN on both sides. To account for the high-temperature operation of the emitter, a temperaturedependent permittivity is used to model the doped Si 14,15 . The optical properties of the aluminum nitride layer are modelled using tabulated permittivity at room temperature 16 .
Using a scattering matrix formalism based on fluctuational electrodynamics, we then calculate the radiative energy transfer between any two layers in the multi-layer structure. In For simplicity, we use the following dimensionless transfer function. Similarly, we denote the transfer function from the active layer of the photovoltaic cell to all the other layers of the photovoltaic cell as → .
Next, to estimate the total energy transfer, QRHT, we calculate the transfer function from the emitter to the whole PV cell denoted as → , where C denotes the whole PV cell. We denote the area of mesa as . This calculation is performed in wavelengths ranging from 100 to 0.7 . Further, to calculate the contributions from the recessed 'rec' region of the emitter, we model it as multi-layer structure as shown in Supplementary Fig. 7b, similar to that followed for the mesa region. The corresponding transfer functions between the rec and active layer, rec and the whole PV cell are denoted by → and → , respectively.

PV cell I-V response modelling
The I-V response of our PV cell is modelled using an equivalent circuit as shown in Supplementary   Fig. 8, where the polarity for the current is defined such that both V and I are positive when the cell generates electricity. We use an experimentally obtained dark I-V response, which is measured under no illumination in vacuum at room temperature, to estimate the parasitic series resistance as 85 based on the slope of I-V curve at large forward voltage and a parasitic shunt resistance as 5 10 Ω.
To calculate the total current generated in our PV cell, we include the current contributions from both the radiative and non-radiative recombination processes. First, we describe the current generated by radiative processes (Irad) as: where e is the elementary charge, TC is the temperature of the cell including the junction, TE is the temperature of the emitter, and subscript denotes all the layers in the PV cell other than the active layer. Here, VFmesa,J and VFrec,J are the view factors associated with both the surfaces, the details of which are discussed in the next section. We assume the internal quantum efficiency of the PV cell to be 100%, or every photon absorbed in the active layer generates an e-h pair. We note that we have used the reciprocal relation for transfer functions, i.e. → → .
Next, to account for current from non-radiative processes, we consider Auger The net recombination current associated with SRH process (ISRH) is given by: where S is the effective surface recombination velocity. The total current (IJ) through the PV cell junction is then given by the sum of all the currents as: .
Further, by accounting for the series resistance and shunt resistance, the total current through the device and the total voltage across the device are: and ⋅ .
The above two characteristic equations are used to obtain the I-V response of the PV cell. The surface recombination velocity 3.5 10 is a fitting parameter adjusted to best match the I-V curve obtained from our model with that of the experimentally measured dark I-V curve.  Supplementary Fig. 9a.

Supplementary
Further, VFrec,J is calculated by approximating the geometry of recessed region and PV cell as two coaxial circular rings with inner diameter 150 μm and outer diameter 190 μm (bottom inset of Supplementary Fig. 9a). An analytical expression for the view factor between circular rings is used as described in ref 18 . The variation of these two factors with gap size is plotted in Supplementary   Fig. 9a. Since our calculation of radiative energy transfer between multilayer structures assumes the structures to be infinitely long in lateral directions, it is important to delineate the effect of view factors in the enhancement of power generation in the near-field. To this end, we plot the variation of theoretically computed PMPP with the gap size. The data presented in Supplementary   Fig. 9b (computed for Temitter = 1270 K) illustrates that inclusion of the corrected view factor results in values of PMPP that are somewhat different from those calculated with a view factor of 1 (infinite plates). In the far-field at 10 μm gap, ~10% variation in the power output could be seen, while the variation is less than 1% at smaller gaps. This is because the energy transfer in near-field is dominated by mesa whose view factor is close to 1 at these gaps as seen in Supplementary Fig. 9a.
Thus, the effect of view factor on power enhancement at smaller gaps is negligible.  Supplementary Fig. 9a, while the green curve corresponds to the assumption of view factor of 1 (infinite plates).

Supplementary Note 9. Characteristics of the PV cell
The current-voltage response of the PV cell under dark conditions (referred to as Dark I-V) is shown in Supplementary Fig. 10a. For this measurement the PV cell was maintained at room temperature in a high vacuum of 1 μTorr. The voltage across the PV cell was varied from -0.5 V to 0.5 V in steps of 10 mV, while the current through the circuit was measured. The variation of fill-factor (FF, Supplementary Fig. 10b) with temperature is plotted along with the theoretical calculations. The small disagreement between measured and calculated FF may likely be due to the uncertainties in our modeling parameters such as the series and shunt resistances.

Supplementary Note 10. Sub-band gap reflection from thin-film PV cell
The improvement in the efficiency of our NFTPV system was achieved by employing a thin-film PV cell as shown in Supplementary Fig. 11a. Typical bulk PV cells suffer from significant subband gap (SBG) absorption as a result of free-carrier absorption in heavily doped layers. The SBG suppression in the thin-film PV cell can be understood by comparing two structures: The current NFTPV system with Temitter = 1270 K at a gap size of 100 nm (Supplementary Fig. 11a) and a thick NFTPV structure, where the n-InP layer was assumed to have a thickness of 200 μm with Temitter = 1270 K at a gap size of 100 nm ( Supplementary Fig. 11b). Using our theoretical model, we then 20 calculate the spectral energy transfers for both cases as shown in Supplementary Fig. 11c as a function of photon energy. The above-band gap flux is similar in both structures because the active layers are of identical thickness. A significant suppression of SBG is seen in a thin-film PV cell as compared to a bulk PV cell, because lower energy photons get absorbed in the 200 μm-thick n-InP layer before getting reflected by the gold layer. Thus, for the thin-film case, low energy photons get reflected back to the emitter improving the overall efficiency. The residual SBG absorption in the thin-film case is due to absorption in the gold layer (imperfect reflector). Figure 11. Thin-film sub-band gap reflection a, Schematic of the thin-film TPV system employed in this study showing the thicknesses of all the layers. b, Schematic of a hypothetical thick TPV structure, where the n-InP layer thickness is increased from 100 nm to 200 μm. c, Spectral energy transfer as a function of photon energy for a thin-film PV cell (solid blue line) and a bulk PV cell (solid red line).