## Abstract

Metamaterial absorbers consisting of metal, metal-dielectric, or dielectric materials have been realized across much of the electromagnetic spectrum and have demonstrated novel properties and applications. However, most absorbers utilize metals and thus are limited in applicability due to their low melting point, high Ohmic loss and high thermal conductivity. Other approaches rely on large dielectric structures and / or a supporting dielectric substrate as a loss mechanism, thereby realizing large absorption volumes. Here we present a terahertz (THz) all dielectric metasurface absorber based on hybrid dielectric waveguide resonances. We tune the metasurface geometry in order to overlap electric and magnetic dipole resonances at the same frequency, thus achieving an experimental absorption of 97.5%. A simulated dielectric metasurface achieves a total absorption coefficient enhancement factor of F* _{T}*=140, with a small absorption volume. Our experimental results are well described by theory and simulations and not limited to the THz range, but may be extended to microwave, infrared and optical frequencies. The concept of an all-dielectric metasurface absorber offers a new route for control of the emission and absorption of electromagnetic radiation from surfaces with potential applications in energy harvesting, imaging, and sensing.

© 2017 Optical Society of America

## 1. Introduction

The study of guided electromagnetic waves began in earnest after the establishment and verification of Maxwell’s equations in the late 1800s. In 1897 Lord Rayleigh investigated the transmission of electromagnetic waves in hollow perfectly conducting tubes [1], while in 1899 Arnold Sommerfeld investigated waves guided by wires [2–4]. In 1909 Sommerfeld’s work was extended by his student - Demetrius Hondros - who obtained a more complete solution to the propagation of surface electromagnetic waves on wires [5]. Hondros went on to work with Peter Debye and they predicted that a solid lossless dielectric of circular cross section could support a TM electromagnetic wave that is confined within the dielectric – except for a small portion lying outside of the cylinder which decreases exponentially [6]. Although experimental verification of the TM guided wave was performed only a few years later [7–9], it would take another two decades before a full mathematical analysis was carried out in 1936 [10]. It was found that both transverse electric (TE) and transverse magnetic (TM) modes exist on a dielectric with a circular cross section, where the electric or magnetic field is transverse to the cylindrical axis. Further, the fields for both the TE and TM modes were shown to be axisymmetric and thus have no azimuthal dependence. Notably, it was also shown that a cylindrical dielectric guide can support hybrid modes which lack azimuthal symmetry, and thus possess angular dependence[10]. As we will show these hybrid modes are of great utility for potential applications, since they are approximate electric dipole (EH) and magnetic dipole (HE) modes, which can be independently tuned through modification of the waveguide geometry [11,12].

Dielectric waveguides have been proposed for many applications including dielectric resonators [13] for bandpass and bandstop filtering applications, antennas, frequency-selective limiters, and various types of oscillators [14, 15]. Dielectric based resonators further have the advantages that they avoid Joule heating – which may occur in a comparable metallic based system – and further may be fashioned from materials which yield good temperature stability for applications. Although traditionally the utility of an all-dielectric waveguide is to transmit electromagnetic waves with little loss, here we form a 2D array of short cylindrical waveguide stacks and use them in the opposite sense. That is, through modification of the geometry, we can tune the lowest order HE and EH modes to overlap in frequency, such that incident energy is not transmitted nor reflected, but rather is completely absorbed entirely within the dielectric. The resulting composite material thus functions as a sub-wavelength metasurface with tunable permittivity and permeability permitting the realization of an all-dielectric absorber [16].

## 2. Design and simulation

We begin by describing the nature of the hybrid modes which are supported by a dielectric cylinder when the k-vector lies parallel to the cylindrical axis (z-axis). A common notation for the two hybrid mode groups is HE and EH, where *H _{z}/E_{z}* « 1 for the HE mode and

*E*« 1 for the EH mode [17]. Thus the EH mode is TE-like since

_{z}/H_{z}*H*is dominant and HE is TM-like since

_{z}*E*is dominant. In order to denote the field variation within the cylinder, three indices are added, i.e.

_{z}*HE*and

_{nml}*EH*. Here

_{nml}*n*indicates the azimuthal variation of the fields and is of the form sin n

*ϕ*and cos n

*ϕ*,

*m*denotes the field variation along the radial direction, and

*l*along the z-axis. Before moving on we note that when the length of the guide is not relevant, and we are only considering the azimuthal and radial modes, we use

*δ*in place of the

*l*-index [17].

We can achieve a highly absorbing state by utilizing the lowest order magnetic dipole (HE_{11}* _{δ}*) and electric dipole (EH

_{11}

*) modes of the cylindrical dielectric waveguide. However, for a single (isolated) dielectric cylindrical waveguide, the HE*

_{δ}_{11}

*and EH*

_{δ}_{11}

*modes do not necessarily possess the same resonant frequency. Notably, the HE mode has zero cutoff and thus may lie at much lower frequencies than the EH mode [11]. However, through understanding of the underlying mechanism responsible for occurrence of the two hybrid modes, we may chose a cylindrical geometry to optimize their frequency overlap. Here we assume that the k-vector of external radiation will be incident parallel to the cylindrical axis – and to the surface normal of a supporting substrate – of our all-dielectric metasurface, and thus only hybrid modes will be excited [18].*

_{δ}We next calculate the minimum height and radius of a dielectric cylinder that supports both electric and magnetic dipole-like hybrid resonant modes at some target frequency *ω*_{0}. Our goal is to achieve a sub-wavelength all dielectric metasurface which realizes high absorption with minimal scattering. We therefore utilize high dielectric materials such that the individual cylindrical resonators are as small as possible – in all dimensions – with respect to the incident wavelength. For a material with a high relative dielectric constant *ϵ _{r}* and for a target frequency with maximum absorptivity of

*A*(

*ω*

_{0}), the wavelength in the guide will be,

The desired minimum height *h* of the guide will support the *l*=1 lowest order mode, and thus can be approximated as a half-wave dipole with,

Although the HE_{11}* _{δ}* mode exists to zero frequency, the EH mode will have a cutoff determined by the radius and dielectric constant of the waveguide. We thus calculate the cutoff radius of a dielectric cylinder that can support the lowest order EH

_{11}

*hybrid mode. We note that although some portion of the wave will lie outside the guide, we want the majority of energy to be contained within the boundary of the cylinder. Thus we constrain the radial k-vector (*

_{δ}*k*) such that the wave undergoes total internal reflection (TIR) inside the cylindrical waveguide. For a guide embedded in vacuum the condition of TIR – in terms of radial wavevector – is then,

_{r}*k*=2

_{g}*π*/

*λ*. The condition placed on

_{g}*k*is used to determine the cutoff radius of our waveguide, which we may find by noting that the solutions to the boundary conditions for modes within the cylinder are determined by roots of the Bessel function of the first kind, i.e.,

_{r}*J*

_{1}(

*k*) = 0, where

_{r}r*r*is the radius of the cylinder [11]. Again we desire a minimum size of our dielectric metasurface, thus we want to solve for the argument of

*J*

_{1}which yields the first nontrivial zero – a value of 3.83. We thus find the cutoff radius is determined by,

For experimental realization of the all-dielectric metasurface absorber, we chose boron doped silicon as the guide material since it realizes a relatively high dielectric value *ϵ _{Si}* = 11.0 with moderate loss at THz frequencies. We choose a target frequency of 1.0 THz and thus, using Eqs. (1)–(4), we calculate the dimensions of our cylindrical waveguide – ignoring the host substrate – to have a height of

*h*= 45

_{Si}*µ*m and a radius of

*r*= 58

_{Si}*µ*m (see Appendix C for more detail). Although we expect an all-dielectric metasurface with the above calculated dimensions to yield high absorption, we first simulate the absorbance for values far away from optimal in order to study the individual hybrid modes and their dependence on geometry.

In Fig. 1(a) we show the absorbance – calculated as A(*ω*)=1-R(*ω*)-T(*ω*) – for a 2D array of cylinders (substrate not included) with a radius of r=40 *µ*m, a height of h=50 *µ*m, and a pitch of p=210 *µ*m. At these dimensions simulation indicates two peaks in the absorption – over the frequency range investigated – at 1.186 THz and 1.294 THz. We find that the lower frequency peak (indicated by the dashed vertical gray line) is due to the magnetic dipole resonance *HE*_{111}, while the higher energy peak (dashed vertical black line) is due to the electric dipole resonance *EH*_{111} [19,20]. In order to investigate the specific dependence of the lowest order HE and EH modes on geometry, we simulate the spectral absorbance as a function of the silicon disk radius from 30 *µ*m to 80 *µ*m, while keeping all other parameters fixed. In Fig. 1(b) we plot a color map of the absorbance as a function of both frequency and radius. We observe that as the radius increases, both the electric (dashed black curve) and magnetic (dashed gray curve) dipole modes shift to lower frequencies. However, simulations reveal that the electric mode is more sensitive to radius which results in a merging of the resonant frequencies as the radius of the disk is increased. We find that near a radius of 60 *µ*m, a highly absorptive mode occurs due to overlap of electric and magnetic resonant frequencies at approximately 1.1THz [16,21]. We also study the dependence of HE and EH modes as a function of height Fig. 1(c) and periodicity Fig. 1(d). The frequency of peak absorbance as a function of height yields a similar dependence as to that of the radius - see Appendix C for more detail. Since our absorbers utilize relatively large dielectric values, the fields fall off rapidly outside of the waveguide. We thus find that *A*(*ω*_{0}) is relatively independent of periodicity and remains large, except for low *p* values where nearest neighbor interactions become significant, and for *p* ≳ *λ*_{0}, where periodic scattering becomes important.

Simulation results presented in Fig. 1 verify the feasibility of achieving high absorption with an array of cylindrical metasurfaces. However we must incorporate a substrate into simulations in order to support our metasurface for eventual fabrication. We chose a Polydimethylsiloxane (PDMS) as a substrate material owing to its low dielectric loss and low thermal conductance. Guided by Eqs. (1)–(4) and simulation results shown in Fig. 1, we select cylindrical dimensions of *h _{Si}* = 50

*µ*m and

*r*= 62.5

_{Si}*µ*m, with a periodicity of

*p*= 172

*µ*m. In Fig. 2(a) we show the simulated spectral dependence of the reflectance (blue curve), transmittance (black curve) and absorbance (red curve). As is evident from Fig. 2(a), just above 1 THz both

*T*and

*R*acquire values near zero and thus simulation indicates that our all-dielectric metasurface achieves a sharp absorption peak at 1.02 THz with a value of A=93.8%. In order to understand the impact of the asymmetry due to the substrate on the hybrid modes, we further investigate the electric and magnetic field distributions at

*ω*

_{0}. In Fig. 2(b) we plot the electric field in the xz-plane, and the magnetic field in the yz-plane at the resonant frequency of 1.02 THz. Fields shown in Fig. 2(b) are plotted at the same phase and we find that they are strongly localized within the dielectric metasurface, with a small portion lying outside (see Appendix A for more details).

## 3. Experimental results

In order to fabricate our dielectric metasurface, a 50-*µ*m thick Boron doped silicon on insulator (SOI) substrate is utilized. A Bosch process was used to pattern the device layer to form an array of silicon disks. The scanning electron microscopy (SEM) image of the patterned SOI sample is shown in Fig. 3(c). The buried 2-*µ*m oxide layer was etched by 49% Hydrofluoric acid to make enough undercut for the following releasing process. Next a 35 *µ*m thick PDMS was spin-coated on a second silicon substrate. The patterned SOI sample was flipped and bonded to the PDMS film by curing at 80 °C for 2 hours. The top SOI handle substrate was then stripped off through application of a lateral force. Finally, the released silicon on the 35-*µ*m PDMS was bonded to a free-standing thick PDMS frame [inset of Fig. 3(c)].

The spectral transmittance and reflectance of the fabricated sample was characterized with a THz time domain spectroscopy (TDS) system. The sample was placed at the focal plane for both reflectance and transmittance setups, and the TDS system was purged with nitrogen gas during the measurement. The absolute reflectance was obtained by referencing measurements with respect to a gold mirror and were made at normal incidence using a beam splitter configuration. Absolute transmittance measurements used an open aperture as a reference and were performed at normal incidence. The measured reflectance (R), transmittance (T) and absorbance (A) are shown in Fig. 3(a) and a peak absorption of A=97.5% at 1.011 THz is realized. In Fig. 3(b) we plot both experimental and simulated *A*(*ω*). Good agreement is evident and the slight difference is due to scattering from the periodic metasurface and surface roughness of the fabricated sample, which is not taken into account in simulation (see Appendix D for more details).

## 4. Discussion

The dipole modes supported by the dielectric cylindrical metasurface extend beyond their surfaces and are thus affected by materials in the surrounding environment. Since a supporting substrate breaks the symmetry, it is not a prior obvious that experimental realization of a dielectric absorber is possible, owing to the necessity of overlapping resonate frequencies of the HE_{111} and EH_{111} modes. In the inset to Fig. 3(b) the spatial distribution of the power loss density is shown in the xz-plane. As can be observed, although the supporting substrate breaks the symmetry, the majority of the power is still dissipated within the silicon metasurface, and only a tiny portion persists in the substrate. However, in order to verify the ability of the all-dielectric metasurface to function as a high absorber of radiation we perform a substrate free simulation and plot A(*ω*) and the power loss density in Fig. 3(d). Simulation reveals that we can achieve a peak absorbance of A=99.6% at a frequency of 1.05 THz.

Our simulated substrate-free metasurface is moderately sub-wavelength and realizes a free space wavelength to particle diameter of *λ*_{0}/2*r*=2.4. The ability of metamaterials and metasurfaces to manipulate electromagnetic waves on a sub-wavelength scales is a salient feature that permits realization of broad-band and wavelength specific absorptivity [22] and emissivity [23]. Like-wise the all-dielectric metasurface absorbers realized here may permit the same high degree of tailorable emission and / or absorption of energy from a surface. Thus in order to quantify the ability of our metasurface to absorb radiation on a sub-wavelength scale, we define an absorption coefficient (*α*) enhancement factor,

*α*and

_{mm}*α*are the absorption coefficients of the metasurface and bulk silicon, respectively, at 1.05 THz (see Appendix B). Our doped silicon yields

_{Si}*α*=42cm

_{Si}^{−1}, whereas our metasurface achieves

*α*=1500cm

_{mm}^{−1}. We thus calculate that the simulated metasurface absorber shown in Fig. 3(d) achieves an absorption coefficient enhancement factor of

*F*=35. We note that for a bulk silicon layer with

_{α}*ϵ*= 11.0 the conventional limit [24] is ${F}_{\alpha}=4{n}_{Si}^{2}=4{\mathit{\u03f5}}_{Si}=44$. However, our dielectric metasurface absorber utilizes significantly less material volume, i.e. the filling fraction is

_{Si}*F*=

*πr*

^{2}

*/p*

^{2}=0.256. Strikingly, our metasurface realizes increased absorption in a reduced volume, and thus we find a total absorption coefficient enhancement factor of [25],

Thus although various dielectric materials may function as absorbers [26,27], they utilize sizes large compared to the wavelength, and / or have a significant portion of absorption occurring within the support substrate.

## 5. Conclusion

We have designed, fabricated and demonstrated an all dielectric THz metasurface absorber. The cylindrical geometry allows independent modification of the EH and HE hybrid modes and we thus realize a highly absorbing state that is localized within the metasurface particles. The absorption mechanism is independent of host substrate and we find a small absorption volume with a total absorption coefficient enhancement factor of 140. Experimentally we realize an absorbance of 97.5% at 1.011 THz and find a good match with simulation. The availability of high thermal insulating properties of dielectric substrates suggest our all dielectric absorber can create a new route in THz thermal imaging/sensing applications. The results demonstrated here are not restricted to the THz frequency range, but like metallic based metamaterials and metasurfaces, may be scaled to operate in nearly any sub-optical band of the electromagnetic spectrum. In contrast to metal-based approaches, an all-dielectric absorber may utilize high melting point and low loss near infrared materials, and thus will be relevant for thermophotovoltaic, imaging, and other sensing applications.

## Appendix A Simulated field distribution at the resonant frequency

As discussed in the main paper, through optimizing the radius and height, both electric and magnetic resonances can be achieved at the same frequency resulting in an impedance matched condition. To further understand how a free-standing silicon metasurface itself can achieve near unity absorption, we simulated the field distributions at the resonant frequency and the results are shown in Fig. 4. According to Fig. 4, at the resonant frequency both electric [Figs. 4(a) and 4(c)] and magnetic [Figs. 4(b) and 4(d)] dipole resonances are excited inside the cylinder and are polarized perpendicular with each other. As is evident from Fig. 4, both the electric and magnetic fields are highly localized within the silicon disk.

## Appendix B Optical properties of bulk silicon and free standing silicon metasurface

In order to verify that the near unity absorption arises from the geometry of the metasuface particles, rather than the bulk material properties itself, we investigated the absorbance of the bare silicon film, while keeping the silicon film thickness the same as the cylinders. As shown in Fig. 5(a), the maximum absorbance of the free standing silicon film is 27.6% at 0.5 THz and decreases to only 14% around 1 THz, which is significantly lower than the near unity absorbance from the cylinders. In addition, we calculated the absorption coefficient of our free standing silicon metasurface (*α _{Meta}*) from the imaginary part of its extracted refractive index (Imag(

*n*)). We plotted the absorption coefficient and the imaginary part of the refractive index of both the bulk silicon and our free-standing silicon metasurface in Fig. 5(b).

_{Meta}## Appendix C Detailed analysis on designing radius and height for perfect absorption with different dielectric constants

In the main paper, we showed that the radius and height of the dielectric cylinder can be determined from the dielectric constant (*ϵ _{r}*) of the cylinder and the designed resonant frequency. According to Eqs. (2)–(4), the height is:

The radius is equal to:

*c*is the speed of light in vacuum and

*f*

_{0}is the target frequency. To show our theoretical calculations generally work for metasurface absorber designs with different dielectric materials, we calculated the radius and height for designing a dielectric metasurface absorber at 1.0 THz with different dielectric constants, and compare with those from the simulations as shown in Fig. 6. In the simulation we kept the periodicity fixed (

*p*= 210

*µm*) and varied the radius and height with the dielectric constant. The calculated results (red solid curves) are close to the simulations (blue open circles). In the calculation, we assume all the fields are within the disk. However, according to simulations shown in Fig. 4, some portion of the fields extend beyond the boundary of the cylindrical particles, and this field leakage results in deviation of the simulated radius from the theoretical calculations. Based on Eq. (8), We fitted the simulations with power law: $a{(x-1)}^{-\frac{1}{2}}$ and the curves are plotted as the red dashed curve in Fig. 6(a). The fitted curve for the radius matches well with the simulated results with the factor

*a*= 200.34

*µm*, and is only 8.7% different from the the calculated factor of 3.83

*c/*(2

*π f*

_{0}) = 182.87

*µm*. Similarly, the height of the cylinder was fit with power law: ${a}^{\prime}{x}^{-\frac{1}{2}}$ and the fitted curve is shown as the red dashed curve in Fig. 6(b). Curve fits match the simulation very well with the fitting factor

*a*′ = 175.58

*µm*, which deviates from the theoretical data of

*c/*(2

*f*

_{0}) = 149.90

*µm*by 14.6%. As a result, the height and radius of the dielectric cylinders from our calculation are close to those from the simulation over a very broad range of the dielectric constants.

## Appendix D Analysis on scattering effect

In our simulations we use unit-cell boundary conditions, where any scattered radiation is collected by Port 1 (reflection) or Port 2 (transmission). Thus, in simulation, scattering would enter into the reflectance (*R _{sim}*) and/or transmittance (

*T*), where

_{sim}*A*= 1 −

_{sim}*R*−

_{sim}*T*. On the other hand, in experiment, the majority of scattered light does not make it back to the detector, and thus the reflectance (

_{sim}*R*) and/or transmittance (

_{exp}*T*) would be smaller. The end result is that

_{exp}*A*= 1 −

_{exp}*R*−

_{exp}*T*is larger than that in simulation. Thus, the scattering (S) may be equal to the difference between the experimental and simulated absorbance, i.e.

_{exp}*S*=

*A*−

_{exp}*A*. For the case studied here, the scattering calculated in this way would be

_{sim}*S*= 97.5% − 93.8% = 3.7%.

We performed a simulation in order to gain more insight into the impact of the scattering. A frequency domain simulation of the exact dielectric absorber structure shown in Fig. 2(a) was performed, and again found the peak absorbance to occur at a frequency of 1.02 THz with a value of A = 93.8%. Our simulation software uses a default emitted port power of 0.5 Watts for each frequency investigated, i.e. at 1.02 THz. We then simulated the distribution of power loss density and integrated it over the dielectric absorber and support dielectric in order to obtain the total absorbed power. We found an integrated power loss of 0.467 Watts to occur within our material and host substrate. The absorbance can be calculated as *A _{power}* = 0.467/0.5 = 93.4%. Thus we find an absorbance with a value that is only 0.4% different from the absorbance evaluated in simulation using

*A*= 1 −

_{sim}*R*−

_{sim}*T*= 93.8%. Although

_{sim}*A*≈

_{power}*A*, it is important to note that they are both unable to determine the amount of scattering. However, in Fig. 3(d) we obtained

_{sim}*A*= 99.6%, for geometrical parameters nearly equal to the case discussed above. Thus since

_{sim}*A*≈

_{power}*A*, the scattering must be negligible since the absorbance determined by power can account for 99.6% of all energy in the simulation, indicating that the scattering is at most 0.4% in this case.

_{sim}We take the above computational analysis to indicate that, if we are able to fabricate a perfect sample, the scattering would be small for the geometrical parameters used for our study. However, in any realistic experimental realization, there would be random errors due to fabricational tolerances and scattering would enter in some amount. The upper bound on scattering in our case is 3.7%.

## Funding

This work was funded by a grant from the Department of Energy (DOE) DE-SC0014372. IVS acknowledges support from the Australian Research Council Discovery Project DP160101585. This work was performed in part at the Duke University Shared Materials Instrumentation Facility (SMIF), a member of the North Carolina Research Triangle Nanotechnology Network (RTNN), which is supported by the National Science Foundation (NSF) (Grant ECCS-1542015) as part of the National Nanotechnology Coordinated Infrastructure (NNCI).

## References and links

**1. **L. Rayleigh, “XVIII. On the passage of electric waves through tubes, or the vibrations of dielectric cylinders,” Phil. Mag. **43**, 125–132 (1897). [CrossRef]

**2. **A. Sommerfeld, “On the propagation of electrodynamic waves along a wire,” Ann. Phys. **303**, 233–290 (1899). [CrossRef]

**3. **K. S. Packard, “The Origin of Waveguides: A Case of Multiple Rediscovery,” IEEE Trans. Microw. Theory Techn **32**, 961–969 (1984). [CrossRef]

**4. **T. K. Sarkar, R. J. Mailloux, A. A. Oliner, M. Salazar-Palma, and D. L. Sengupta, *History of Wireless* (John Wiley & Sons, 2006). [CrossRef]

**5. **D. Hondros, “On electromagnetic wire waves,” Ann. Phys. **335**, 905–950 (1909). [CrossRef]

**6. **D. Hondros and P. Debye, “Electromagnetic waves in the electric wire,” Ann. Phys. **337**, 465–476 (1910). [CrossRef]

**7. **H. Zahn, “On the detection of electromagnetic waves on dielectric wires,” Ann. Phys. **354**, 907–933 (1916). [CrossRef]

**8. **H. Ruter and O. Schriever, “Electromagnetic waves on dielectric wires,” Schriften des Naturalwissenschaftlichen vereines fur Schleswig-Holstein **16**, 2 (1916).

**9. **O. Schriever, “Electromagnetic waves on dielectric wires,” Ann. Phys. **368**, 645–673 (1920). [CrossRef]

**10. **J. R. Carson, S. P. Mead, and S. A. Schelkunoff, “Hyper - Frequency Wave Guides - Mathematical Theory,” Bell Syst. Tech. J. **15**, 310–333 (1936). [CrossRef]

**11. **E. Snitzer, “Cylindrical dielectric wavegnide modes,” J. Opt. Soc. Amer. **51**, 491–498 (1961). [CrossRef]

**12. **Y. Kobayashi and S. Tanaka, “Resonant Modes of a Dielectric Rod Resonator Short-Circuited at Both Ends by Parallel Conducting Plates,” IEEE Trans. Microw. Theory Techn. **28**, 1077–1085 (1980). [CrossRef]

**13. **R. D. Richtmyer, “Dielectric resonators,” J. Appl. Phys. **10**, 391–398 (1939). [CrossRef]

**14. **J. K. Plourde and C. L. Ren, “Application of Dielectric Resonators in Microwave Components,” IEEE Trans. Microw. Theory Techn. **29**, 754–770 (1981). [CrossRef]

**15. **P. G. Darko Kajfez, *Dielectric Resonators* (SciTech Publishing, 1998), 2.

**16. **M. A. Cole, D. A. Powell, and I. V. Shadrivov, “Strong terahertz absorption in all-dielectric Huygens’ metasurfaces,” Nanotechnology **27**, 424003 (2016). [CrossRef]

**17. **R. K. Mongia and P. Bhartia, “Dielectric resonator antennas - a review and general design relations for resonant frequency and bandwidth,” Int. J. Microwave Mill. **4**, 230–247 (1994).

**18. **C. Yeh and F. Shimabukuro, *The Essence of Dielectric Waveguides* (Springer, 2008). [CrossRef]

**19. **J. V. D. Groep and a. Polman, “Designing dielectric resonators on substrates : Combining magnetic and electric resonances,” Opt. Express **21**, 1253–1257 (2013).

**20. **Q. Zhao, J. Zhou, F. Zhang, and D. Lippens, “Mie resonance-based dielectric metamaterials,” Mater. Today **12**, 60–69 (2009). [CrossRef]

**21. **I. Staude, A. E. Miroshnichenko, M. Decker, N. T. Fofang, S. Liu, E. Gonzales, J. Dominguez, T. S. Luk, D. N. Neshev, I. Brener, and Y. Kivshar, “Tailoring directional scattering through magnetic and electric resonances in subwavelength silicon nanodisks,” ACS Nano **7**, 7824–7832 (2013). [CrossRef] [PubMed]

**22. **X. Liu, T. Starr, A. F. Starr, and W. J. Padilla, “Infrared spatial and frequency selective metamaterial with near-unity absorbance,” Phys. Rev. Lett. **104**, 207403 (2010). [CrossRef] [PubMed]

**23. **X. Liu, T. Tyler, T. Starr, A. F. Starr, N. M. Jokerst, and W. J. Padilla, “Taming the blackbody with infrared metamaterials as selective thermal emitters,” Phys. Rev. Lett. **107**, 045901 (2011). [CrossRef] [PubMed]

**24. **E. Yablonovitch, “Statistical ray optics,” J. Opt. Soc. Am. **72**, 899–907 (1982). [CrossRef]

**25. **Z. Yu, A. Raman, and S. Fan, “Fundamental limit of nanophotonic light trapping in solar cells,” Proc. Natl. Acad. Sci. U.S.A. **107**, 17491–17496 (2010). [CrossRef] [PubMed]

**26. **Y. Z. Cheng, W. Withayachumnankul, A. Upadhyay, D. Headland, Y. Nie, R. Z. Gong, M. Bhaskaran, S. Sriram, and D. Abbott, “Ultrabroadband Plasmonic Absorber for Terahertz Waves,” Adv. Opt. Mater. **3**, 376–380 (2015). [CrossRef]

**27. **Y. Sheng, J. Zhu, W. Xu, W. Jiang, J. Yuan, G. Yin, L. XIe, Y. YIng, and Y. Ma, “High-performance terahertz wave absorbers made of silicon-based metamaterials,” Appl. Phys. Lett. **107**, 073903 (2015). [CrossRef]