Imaging Metasurfaces based on Graphene-Loaded Slot Antennas

Spectral imagers, the classic example being the color camera, are ubiquitous in everyday life. However, most such imagers rely on filter arrays that absorb light outside each spectral channel, yielding ~1/N efficiency for an N-channel imager. This is especially undesirable in thermal infrared (IR) wavelengths, where sensor detectivities are low, as well as in highly compact systems with small entrance pupils. Diffractive optics or interferometers can enable efficient spectral imagers, but such systems are too bulky for certain applications. We propose an efficient and compact thermal infrared spectral imager comprising a metasurface composed of sub-wavelength-spaced, differently-tuned slot antennas coupled to photosensitive elements. Here, we demonstrate this idea using graphene, which features a photoresponse up to thermal IR wavelengths. The combined antenna resonances yield broadband absorption in the graphene exceeding the 1/N efficiency limit. We establish a circuit model for the antennas' optical properties and demonstrate consistency with full-wave simulations. We also theoretically demonstrate broadband ~36% free space-to-graphene coupling efficiency for a six-spectral-channel metasurface. This research paves the way towards compact CMOS-integrable thermal IR spectral imagers.

We take spectral imaging for granted in daily life. Our eyes are spectral imagers, providing information about the composition of what we see. The infrared electromagnetic band covers many chemical absorption resonances and thus also reveals compositional information. In particular, infrared spectral imaging is applied in areas such as gas emission monitoring, 1-7 ecological monitoring, 8-10 food quality control, [11][12][13] waste sorting, 14 biological research 15 and oceanography. 16 Spectral imaging aims to measure a "data cube" representing light intensity over two spatial dimensions x and y and one spectral dimension λ with N channels. Scanning spectral imagers sequentially measure different portions of the data cube over multiple exposures to form the full data cube. A common example is the pushbroom scanner which measures x × λ data cube slices while scanning y and is thus typically associated with satellites and conveyor belts where either the camera or subject is gradually moving in one direction. 1,12,14,17 The spectral axis may also be scanned such as in tunable filter-based imagers, 3,18 which feature at most 1/N light utilization efficiency, or Fourier transform interferometer-based imagers which are bulky and have moving parts. 2 In contrast, snapshot spectral imagers (SSIs) capture a data cube with a single exposure. 19 This may be achieved using a color filter array similar to that of a color camera, thus limiting efficiency to 1/N . [19][20][21] Another category of SSIs uses dichroic or dispersive optics to break up incoming light by wavelength before arriving on a focal plane array (FPA). There are many variations of this approach, 15,[22][23][24][25] but they all require increasing the etendue of the incoming light beam by a factor of N , leading to an unfavorable tradeoff between total FPA area, input acceptance angle and spectral resolution. 19 Compared to these technologies, the category of imagers based on inherently multispectral pixels is less explored. One such example is the Foveon RGB sensor, which extracts three different electrical signals from different depths in the optically active silicon region, as shorter wavelengths are absorbed closer to the sensor surface. 26 Another approach uses nanoantennas, optical resonators of subwavelength dimensions which nevertheless feature absorption cross sections of order λ 2 if the antenna is conjugate impedance matched with its load; or, equivalently, if the antenna is critically coupled to the vacuum. Here, we propose a thermal IR multispectral imager where N differently sized metallic slot antennas with infrared-sensitive loads targeting N spectral channels are tiled to form a metasurface featuring efficient free space-to-load optical energy transfer. Figure 1a shows such a metasurface for N = 6. We model graphene as the photosensitive load because its broadband absorption in the mid-IR 27 and processing flexibility 28,29 make it suitable for this platform. Not only do the antennas sort incident light by spectral channel, but they also enhance the absorption of the graphene load, bridging the gap between the impedance of free space and graphene, which, when undoped, has an optical sheet resistance no lower than 16.1 kΩ. 27 Figure 1b shows in detail a single such antenna-coupled graphene photodetector. This detector is designed for a strong photothermoelectric response, in which absorbed light heats up the electron gas in the graphene, resulting in an electromotive force due to the Seebeck effect. The graphene channel is assumed to be isolated from the metasurface by a severalnanometer layer of dielectric, thin enough to not impact the optical properties of the system. The asymmetric position of the graphene channel with respect to the slot allows half of the graphene channel to be gated by metal underneath, yielding the asymmetric graphene Fermi level profile needed for a nonzero net photoresponse. 30,31 Note that while perfect absorption in this wavelength range has been demonstrated in heavily doped graphene accompanied with metal nanostructures, 32 we limit our consideration to undoped graphene as the peak Seebeck coefficient occurs at very low doping levels. 30 For optical absorption, slot antennas offer a few advantages over planar designs such as dipole or bowtie antennas. They have unidirectional radiation patterns, and thus an array of them can perfectly absorb an incident beam, whereas planar antennas require a quarter-wave back-reflector to do so. 33,34 The wavelength dependence of the back-reflector phase complicates design of broadband absorbing metasurfaces and exacerbates undesirable antenna-antenna coupling. Additionally, since planar antennas must be supported by transparent dielectric, they cannot be embedded in a CMOS process as the inter-layer dielectric strongly absorbs thermal infrared radiation, 35,36 whereas for slot antennas the dielectric on top and in the perforations may be etched away without sacrificing the mechanical integrity of the antenna.

Results and Discussion
Z gr Z A Figure 2: a) Graphene-loaded slot antenna with physical features corresponding to the components in circuit b) labelled. b) Circuit schematic of a graphene-coupled slot antenna. V A represents incoming light, Z A is the radiation impedance of the slot aperture, Zgr is the impedance of the graphene sheet and Zwg is that of the slot, effectively a short-circuit waveguide stub.
We model the slot antenna depicted in Figure 2a as an aperture antenna fed by a rectangular waveguide terminated in a short circuit a distance d behind the aperture. We represent the graphene sheet as a shunt impedance connected in parallel with the waveguide stub. Figure 2b illustrates this circuit with a Thévenin equivalent radiation impedance ZA and source VA representing the aperture antenna. 37 The rectangular waveguide stub presents an impedance where Z0 and n eff are the characteristic impedance and effective index of the TE10 mode of the slot and k0 is the vacuum wavenumber. r is the Fresnel reflection coefficient between vacuum and metal for an s-polarized plane wave at an incident angle of arccos(n eff ), which describes the TE10 mode.
The graphene presents a mostly resistive impedance of using power/current impedance normalization. 38 w and h are defined in Figure 1b and σgr is the sheet conductance of the graphene, modeled here as intrinsic. We calculate ZA using finite element simulations for various h and w; we provide more details on these calculations in the Methods section. See Supplementary Figure 1 for an example of the frequency dependence of the impedances in the circuit. Define ηgr as the fraction of available power from Thévenin source VA, ZA that is dissipated in Zgr. Solving the above circuit, we arrive at where represents reciprocal addition. Define Agr = Pgr I inc as the partial absorption cross section of light of intensity Iinc coupled into the graphene and Pgr as the power absorbed in the graphene. For a lossless, conjugate matched antenna, antenna theory predicts Agr,max = D(θ,φ) 4π λ 2 where D (θ, φ) is the antenna's directivity at the given incident angle and polarization, 37 which we calculate using finite element simulations. θ here is the polar angle and φ is the azimuthal angle. We omit the polarization angle in our notation, implicitly setting it to maximize D. Additional ohmic losses and impedance mismatch then reduce the actual absorption cross-section into graphene by a factor ηgr, i.e.
We obtain Agr from FDTD simulations of plane waves incident on the graphene-loaded antennas, which we then use to calculate ηgr via Equation 4. Fig. 3 compares ηgr between the model described by Equation 3 and FDTD results for antennas of various dimensions. The data show that the model is accurate to within 10% of the ηgr peak amplitude and 2% of the resonance wavenumber. We attribute these discrepancies to aspects not captured by the quasi-analytical model, such as our assumption of a perfectly conducting outer antenna face and finite meshing. Despite these shortcomings, the present model allows us to predict slot antenna absorption properties to tolerances comparable to the uncertainty due to variations in metal quality.
To further validate our model, we artificially scale the sheet conductance of the graphene load by factors ranging from 0.33 to 5 and compare the resulting ηgr between the model and FDTD for a single such antenna. The results,  shown in Figure 4, show that our model accurately predicts the sublinear scaling of ηgr with respect to the load conductivity, with the ηgr peak amplitude reaching 0.8 for a load conductivity of 5σgr.
Having modeled the individual components, we now discuss broadband absorbing metasurfaces incorporating differently tuned antennas tiled in a periodic array. We use a three-step process to design such metasurfaces. We first constrain d and w to be the same for all antennas in the metasurface, and choose their values to yield high peak ηgr for antennas resonant in the targeted wavelength range. Secondly, we choose the values of h for the antennas, following the heuristic that at the wavelength where one antenna's ηgr falls to half its peak amplitude, the next antenna's ηgr should have risen to half its peak amplitude. Finally, we choose the arrangement and pitch of the antennas to be as closely packed as possible while satisfying qualitative fabricability constraints. We also avoid juxtaposing antennas of adjacent wavelength channels to minimize antenna-antenna crosstalk.
The antenna pitch, more accurately described by the Bravais lattice vectors of the array, is a critical parameter in determining the potential absorption efficiency of the array. Light incident from a given direction can only be scattered by the two-dimensional diffraction orders of the lattice. The array can only perfectly absorb an incoming light beam if no nonzero diffraction orders fall within the light cone, barring the event that all higher diffraction orders overlap with nodes in the individual unit cell radiation pattern. By "light cone", we refer here to the region in the Fourier transform space of the xy-plane for which radiation can occur, namely k 2 x + k 2 y < k 2 0 . For a square lattice, if the lattice pitch a < λmin/2, no higher diffraction orders are within the light cone for any incident angle. In practice, the limited numerical aperture of imaging systems relaxes this constraint.  Targeting the 6 − 10 µm wavelength range, we follow the above methodology and arrive at an array of six antennas with d = 2.25 µm, w = 400 nm, and h = 3.41, 3.81, 4.41, 5.04, 5.80, 7.36 µm, uniformly spaced and tiled as shown in Figure 1a with a 7 µm by 12.667 µm unit cell. Figure 5 shows the absorption efficiency of normally incident light into the graphene load of each antenna as well as their sum, simulated with FDTD. With an average efficiency of 36% across the 1050 cm −1 to 1600 cm −1 band, this structure improves upon the 1/N limit of filter array-based imagers by roughly a factor of two. Note that the unit cell highlighted in Figure 1a is not the primitive unit cell of the lattice, although it was used as the FDTD simulation region due to software constraints.
To further understand the physics of these metasurfaces, we vary the pitch of the unit cell while keeping all other parameters constant. The resulting total absorption efficiency spectra for six different cases are shown in Figure 6a efficiency decreases with increasing unit cell size. This can be understood by analyzing the diffraction characteristics of the various metasurfaces. Sparser metasurfaces have tighter reciprocal lattices and thus more diffraction orders are available within the light cone. For normally incident illumination as we are using here, the minimum wavelength for which no higher-order diffraction occurs (i.e., completely specular reflection) is given by λc = max a −2 , where a1 and a2 are the unit cell dimensions. The corresponding wavenumbers for the various unit cells used here are listed in Figure 6c. For the 9 µm and 11 µm width unit cells, the threshold wavenumber falls in the middle of the range of interest, resulting in reduced efficiency as diffracted light cannot participate in destructive interference with specularly scattered light. This is especially apparent for the 9 µm unit cell, which exhibits a clear spectral transition between high and low efficiency at the diffraction threshold. Note that choosing an excessively tight antenna spacing is also detrimental, not only due to fabrication difficulty, but also because graphene detectors feature minimum Johnson noise-dominated noise-equivalent power for channel lengths comparable to the hot carrier cooling length, which can range from 100 nm to over 1 µm depending on the substrate and graphene quality. 30,31,39 Exceedingly tight antenna spacings may not provide room for such long graphene channels.
For a spectrally sensitive metasurface to be practical, not only must it maintain a high absorption efficiency for a reasonable range of incoming light directions, but also the absorption spectra of the individual antennas must not shift or distort too strongly as the incoming light angle varies. The directional dependence of our metasurface arises from two factors: The directivity profile D (θ, φ) of the individual antennas, and array effects resulting from interference and antenna-antenna coupling. Although full angle-dependent simulation results for our gold metasurfaces are outside the scope of this paper due to the extreme computational overhead of off-angle periodic structure simulations, 40 we can still provide insight by elaborating on the aformentioned factors, and we also perform off-angle simulations of a simplified metasurface constructed of Perfect Electrical Conductor (PEC) which permits a much larger mesh size. Supplementary Figures 3a, b and c display the directivity profile of a 6.5 µm × 400 nm antenna on resonance. This profile is similar to those of the other antenna lengths used in the metasurface. Intuitively, the directivity decreases as the incident angle approaches the long axis of the antenna, and we thus expect a similar trend in the directional dependence of the array. In Supplementary Figure 3d, we plot for three wavenumbers the sets of incident light directions, represented as components kx, ky of the incident wavevector k, for which only specular reflection from the metasurface occurs. For λ −1 = 1050 cm −1 , all light incident within 45 • of normal is only specularly reflected. This range decreases with increasing wavenumber until normally incident light is pinched off at 1579 cm −1 . As in Figure 6, we expect efficiency to suffer when the specular reflection-only condition is not met. Supplementary Figure 4 shows the results of offangle simulations of the simplified PEC metasurface. The data show that for light incident off-angle but perpendicular to the antennas' long axes, antenna resonances falling outside the specular reflection-only region are subject to decreased absorption efficiency as well as blue-shifting. For light incident off-angle and perpendicular to the antennas' short axes, we obtain similar results, except that the peaks red-shift instead of blue-shift, and we observe an overall reduction of the absorption at steep incident angles due to the reduced directivity.
We also investigate metasurfaces comprising numbers of spectral channels besides N = 6. Supplementary Figure 5 shows the geometric details and simulation results for metasurfaces with N = 3, 4, 5, and 8 as well as the default value of 6. We achieve good results with uniformly high absorption efficiency for N = 5 and 6. For lower values of N , the wide frequency spacing between the individual resonances yields deep troughs in the overall absorption efficiency curve, whereas for higher values of N , excessive overlap between the antenna resonances causes the overall metasurface efficiency to suffer.
With realistic metal, the efficiency of these devices is ultimately limited by ohmic losses. However, more advanced antenna designs can be used to improve ηgr. As it turns out, not only can the copper layers in the back end of a CMOS process be designed to incorporate slot antennas, but they also provide a convenient medium to realize those antenna designs that would be exceedingly difficult to fabricate in an academic setting. 43 Figure 7a introduces such a design in which the slot inlet is narrower than the internal cavity width. This design concentrates the electric field around the graphene, which reduces Zgr without increasing the TE10 mode loss. Figure 7b shows a CMOS adaptation of this de-  Figure 1b), while "slitted" refers to the slitted design in a). "Evap" refers to evaporated gold as used throughout the paper, 41 whereas "SC" refers to single crystal gold. 42  sign, in which the walls of the cavity are perforated to comply with via design rules. If the perforations are sub-cutoff for the resonant wavelength and sufficiently deep, they do not leak light. Additionally, it would be necessary to etch away the inter-layer dielectric to prevent light absorption.
Besides different antenna designs, material quality also affects efficiency considerably. We explore both of these variations in Figure 7c, which plots ηgr versus wavenumber for intrinsic graphene as simulated by FDTD for narrowed inlet ("slitted") antennas and normal slot antennas, as well as with a single-crystal gold model 42 in addition to our default evaporated gold model. The data show that adopting a slitted antenna design increases the peak ηgr from 0.4 to 0.6, and then to 0.9 for single-crystal gold. We can thus hope to achieve ηgr ≈ 0.6 for CMOS-integrated antennas, as copper has been shown to exhibit slightly superior plasmonic properties to gold given suitable deposition conditions. 44 While one can achieve high efficiencies with clever antenna designs and more opaque loads than monolayer graphene, the Q is ultimately bounded by that of a sealed metal cavity which we simulate to be about 40 for the present gold model and antenna shapes. Silver has less mid-IR optical loss than copper or gold, but unlike copper it is not considered CMOScompatible and thus may be difficult to integrate. Polar materials supporting optical phonons in the mid-IR are reported to have high Q plasmonic resonances and are thus worth investigating if higher Q is necessary, although plasmonic behavior only occurs in narrow wavelength bands, limiting applicability. 45 We can apply data describing the resistivity of gated graphene to our model to estimate the room temperature detectivity of the spectral imager described by Figures 1 and 5. Using the methods described in Song et al. 31 and the electrical properties of polycrystalline, non-annealed graphene achieved by de Fazio et al., 46 we arrive at the values in Table 1 with and without accounting for graphene contact resistance for normally incident x-polarized light. Using the more advanced antenna designs shown in Figure 7, the detectivities would reach the 10 8 Jones range. We discuss these calculations in the Methods section. For comparison, more conventional bolometer-based thermal IR FPAs operating at room temperature have been reported to achieve detectivities in the 1-2 × 10 9 Jones range, although such devices do not feature spectral resolution and are limited to millisecond-range response times. 47,48 Finally, we would like to emphasize the general applicability of the antenna metasurface concept to other wavelength ranges, photosensitive elements and antenna designs. Indeed, Tamang et al. have proposed a similar concept applied to RGB imaging where silicon nanorods act as both the antenna and sensitive element. 49 We propose that besides graphene and other 2D materials, III-V or HgCdTe semiconductor photosensitive elements could also be incorporated by a transfer printing heterointegration process. 50 The slot antenna-based metasurface imager approach could also scale to terahertz, where Ohmic losses are much less than in the mid-IR 51 and the antennas could be fabricated directly in a circuit board-like platform.

Conclusion
In conclusion, we introduced a six-spectral-channel graphene-coupled slot antenna metasurface with 36% efficiency functioning as a spectral imager in the thermal IR, as well as a model for estimating the optical properties of the individual antennas therein. This device is appropriate for integration in the wiring layers of a CMOS process with suitable post-processing to remove inter-layer dielectric within the cavity and transfer graphene. We have shown that more sophisticated antenna designs can improve the efficiency of optical energy transfer to the load to above ηgr = 0.6. Further research on this concept may focus on experimental demonstration of the absorption enhancement functionality, or on optimizing the device design to meet certain engineering goals.

Simulation details
Unless otherwise specified, we use the evaporated gold optical model described in Palik et al. throughout the paper. 41 We model graphene as an infinitely thin conductive sheet using the optical conductivity model described in Hanson. 52 As input parameters to the model, we use a temperature of 300 K, intrinsic graphene (zero Fermi level), and a scattering rate Γ = 0.514 meV.
We use the Lumerical FDTD package for our FDTD simulations. For simulations of individual antennas, we use x-, y-and z-meshes of 8 nm in the vicinities of the slot aperture and slot bottom, as well as a z-mesh of 40 nm within the slot. As such, the finest meshes coincide with metallic surfaces and corners, allowing us to capture the nonzero skin depth of the metal, whereas the more gradual z-dependence of the fields inside the slot permits a coarser mesh. We find that a minimum mesh size of 8 nm yields converged results for these simulations. We use PML boundary conditions on all sides except for the −z side, where we use a metallic boundary condition as the light does not penetrate beyond the slots anyway. We also apply symmetry conditions across the xz and yz planes. For the metasurface simulations, we use the same meshing scheme, but with a fine mesh of 15 nm and a z mesh of 100 nm within the slot due to computational resource availability limits. To illustrate the error associated with this choice of mesh, we plot mesh-dependent absorption efficiency curves for a 9 µm by 13.333 µm unit cell, N = 6 metasurface in Supplementary Figure 6, with the mean efficiency averaged between 1050 cm −1 and 1600 cm −1 plotted in the inset. The results validate our qualitative conclusions and put an approximately 3% relative error bound on the spectrally averaged efficiencies of the metasurfaces, although the coarse mesh does somewhat distort the actual spectra. For the metasurface simulations, we change the x-and yboundary conditions to Bloch boundary conditions to reflect the periodic nature of the metasurface, and we apply symmetry across only the xz plane.

Impedance model details
We calculate Z0 and n eff for our impedance model using the Lumerical MODE waveguide mode solver with an 8 nm mesh. We use the graphene model from Falkovsky, 27 which gives almost identical results to the Hanson model used by Lumerical for the parameters we use. We use Ansys HFSS finite element software to calculate ZA. These simulations excite the aperture from within by its TE10 mode yielding the S11 scattering coefficient of the internally reflected light, from which the software calculates ZA. In the finite element simulations, we model the aperture and slot as perfect electrical conductors, as we expect the real part of the antenna impedance to be dominated by radiative loss (rather than ohmic loss) and the imaginary part by energy storage in the near field of the aperture (rather than plasmonically within the metal). From these same simulations we also extract the antenna directivity D (θ, φ).

Detectivity estimation
We base our estimation of the detectivity D * on the formulation put forth in Song et al. 31 To calculate the electronic temperature profile of graphene suspended over the slot, we solve the 2-dimensional partial differential equation: Here κ represents the electronic planar thermal conductivity of the graphene; ∆T el is the difference between the thermally excited electronic temperature T el and the lattice temperature T0 = 300 K, γ represents the electron-phonon thermal decay rate, and C el represents the electronic heat capacity of graphene. α is the efficiency with which optical energy absorbed by the graphene is deposited into the electronic system on a sub-picosecond timescale, taken to be unity since the incident photon energy is below graphene's optical phonon energy. 0Ṅ represents the intensity profile of absorbed light. j is the electrical current density, and Π refers to the Peltier coefficient. We choose an antenna with w = 400 nm and h = 5.5 µm in these calculations, extracting the spatial dependence of 0Ṅ from Ansys HFSS finite element simulations. For the graphene's conductivity σ as a function of Fermi level, we use the data measured by de Fazio et al. for unannealed, polycrystalline graphene; 46 this is then used to calculate κ via the Wiedemann-Franz law and Π as well as the Seebeck coefficient S via the Mott formula. The value of γC el is estimated by assuming a electronic thermal cooling length of κ/γC el = 1 µm, an empirical value. 39 As shown in Figure 1b, the graphene is assumed to be terminated at the slot edge on one side, and is taken to extend 400nm past the slot edge on the other side, where its Fermi level is gated through the metal to the n-type Seebeck coefficient peak. The graphene Fermi level in the suspended region is simply taken to be the zero-gate-voltage Fermi level from de Fazio et al. as it cannot be controlled. For simplicity we assume a sharp jump in the values of σ, κ, Π and S between the n-and p-doped sides of the device, neglecting fringing fields from the gate. The graphene channel is assumed to be short-circuited with graphene-metal contact resistances Rc of either 0 or 1000 Ω µm per contact, the latter being consistent with 1-dimensional contacts to graphene near the Dirac point. 53 The average ∆T el at the graphene p-n junction and the Seebeck coefficients on either side determine the thermal electromotive force via E = −S∇T , 54 which in turn determines j via the total device resistance. Thus, Equation 5 including the Peltier term may be solved directly, as j is a linear functional of ∆T el . Solving for ∆T el over two spatial dimensions, we find linear thermal decay profiles along the +x and −x directions away from the ∆T el peak which indicates that the device is in the short-channel regime where carrier cooling is dominated by the ∆T el = 0 boundary conditions, and κ/γC el is large enough for the electron-phonon interaction term to be inconsequential.
Having obtained the short-circuit responsivity R under zero bias as such, we calculate the noise-equivalent power of the device assuming Johnson noise at 300 K as the dominant noise source, a reasonable assumption for an unbiased device. 55 The detectivity for the array is calculated incorporating the antenna pitch, noting that there are two antennas per unit cell. To account for the decreased optical absorption efficiency of a metasurface loaded with graphene doped to the Seebeck coefficient peaks of roughly ±0.05 eV, we redo the simulation used to generate Figure 5 with the graphene doped as such. We plot the resulting absorption spectra in Supplementary Figure 7, which shows a mean absorption efficiency of 33% averaged between 1050 cm −1 and 1600 cm −1 . We finally calculate the external values of responsitivity, NEP, and detectivity by scaling the internal values by this efficiency factor.

Supporting Information
Additional plots to supplement the main text, including:

Radiation angle θ [degrees]
Supplementary Figure 3: Different representations of the on-resonance far-field radiation pattern of a 6.5 µm-long slot antenna, which is qualitatively very similar to those of antennas of other lengths. a) 3D-representation of the far-field directivity pattern. D (θ, φ) goes to zero for light incident along the long axis of the aperture. b) Directivity vs. elevation angle θ along the azimuth of the antenna's long axis, φ = 90 • . c) Directivity in directional cosine space. Note also that k x = k 0 u and k y = k 0 v where k x and k y are the x-and y-components of the incident wavevector. d) Depiction in k x k y -space of the light cone edges (colored rings) and the specular reflection-only incident light directions (colored patches) for three different wavenumbers for the 7 µm × 12.667 µm 6-antenna unit cell design in Main Figure 1a. The black dots represent the reciprocal lattice points of the metasurface. For each wavelength, the specular reflection-only region is the set of remaining {k x , k y } after subtracting from the light cone all copies of the light cone transposed by all nonzero reciprocal lattice vectors.    Figure 5, but with the graphene doped to 0.05 eV which decreases the graphene sheet conductivity to a degree, especially for longer wavelengths. Here, the mean efficiency averaged between 1050 cm −1 and 1600 cm −1 is 33%. This doping level is chosen to maximize the Seebeck coefficient for the graphene model used to approximate the device detectivity as discussed in the main text.