## Abstract

Arrays of dielectric cylinders support two fundamental dipole active eigenmodes, which can be manipulated to elicit a variety of electromagnetic responses in all-dielectric metamaterials. Dissipation is a critical parameter in determining functionality; the present work varies material loss to explore the rich electromagnetic response of this class of metasurface. Four experimental cases are investigated which span electromagnetic response ranging from Huygens surfaces with transmissivity T = 94%, and phase *ϕ _{S}*

_{21}= 235°, to metasurfaces which absorb 99.96% of incident energy. We find perfect absorption to be analogous to the driven damped harmonic oscillator, with critical damping occurring at resonance. With high phase contrast, transmission, and absorption all accessible from a single system, we present a uniquely diverse all-dielectric system.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Metamaterials were first demonstrated in a – now common – metal and insulator based paradigm. This ubiquitous design has contributed to a burgeoning field in which a variety of tailored electromagnetic responses have been achieved from negative refraction, to cloaking, and perfect absorption over a broad spectrum [1–6]. Many aspects of the metal based metamaterial architecture are of continued interest – both fundamentally and in application – however, the convention has inherent limitations. Technologies with high functional temperatures (>1000°*C*), such as thermophotovoltaics (TPVs), operate above the melting point of most common metals [7], and thus cannot be implemented with metals. To this end, metamaterials consisting entirely of dielectric materials present a solution, and with their relatively recent emergence, [8–10] the full potential is yet to be realized.

All-dielectric metamaterials (DMMs) have seen an increased popularity in recent years. Several studies have demonstrated resonant transmission, reflection, and absorption, showcasing the versatility of the dielectric based system. Active tuning is also possible, and approaches ranging from mechanical deformation to use of liquid crystals have been demonstrated [11,12]. Investigations into the underlying physics of DMMs have advanced characterizations varying from Mie theory [13], to multi-pole expansions [14], and waveguide theory [10]. An initial grasp of the mechanisms, and some design constraints, have emerged for resonant particle arrays, however a complete description is absent, and critical portions of the parameter space remain unexplored. Prior studies of periodic metallic metamaterials and surface phonon-polariton gratings have shown that absorbing systems depend strongly on a balance of radiative loss rates with dissipative rates [15–17]. From these and similar studies on periodic absorbing systems, two phenomena can be identified: dielectric loss is linked to the dissipative rate, and filling fraction to the radiative rate [15], while radius and height set the mode frequencies in the system [10]. The resonators’ material properties, and filling fraction will be shown to be critical parameters that govern the electromagnetic properties, once a system has been chosen, and can be thought of as parameters of dissipation and radiation, respectively.

The current work is focused on loss in the dielectric metasurface (DMS), and offers insight into the versatility of the all dielectric metamaterial. We experimental demonstrate the impact of material loss on the function of the DMS, and a broader exploration of the full parameter space is performed via simulation. The effects of loss on the scattering properties of the DMS are investigated utilizing temporal coupled mode theory coupled with effective material parameters obtained by numerical simulations. Temporal coupled mode theory allows the system to be described by a radiative and dissipative rate, the ratio of which we show to be proportional to dopant concentration in a silicon material system. The special case of perfect absorption is found to parallel the damped harmonic oscillator, and is intuitively understood within the context of the effective optical constants. Results suggest that modification of loss in the DMS is able to elicit a wide ranging response without any change in geometry, with a single DMS achieving a diverse set of fundamental electromagnetic properties.

## 2. Huygens metasurface design

The all-dielectric Huygens surface has been demonstrated in a number of geometries [9]. Most systems studied to date use simulation to achieve appropriate geometries and subsequently attribute the response to Mie resonance. The present work focuses on cylindrical resonators in a square periodic array, following recent work on dielectric metamaterial absorbers [10, 18]. In the work by Liu *et al.* [10] the metasurface is found to achieve perfect absorption when two waveguide modes (*HE*_{111}, *EH*_{111}) are supported simultaneously in the unit cell. Both modes are fundamental in the cylindrical waveguide and cutoff conditions are derived as,

*r*is radius,

*J*

_{1,1}is the first positive zero of the Bessel function of the first kind,

*λ*is the guided wavelength given by ${\lambda}_{g}={\lambda}_{0}\sqrt{{\mathrm{\u03f5}}_{1}}$. Here

_{g}*ϵ*

_{1}is the real part of the complex permittivity $(\tilde{\mathrm{\u03f5}}={\mathrm{\u03f5}}_{1}+i{\mathrm{\u03f5}}_{2})$ of the cylinder,

*h*is the height of the cylinder, and

*λ*

_{0}is the wavelength in free space at resonance. Although, the work by Liu

*et al.*establishes an analytic guide for design of the all dielectric metasurface absorber, we show here that design by Eqs. (1) and (2) is in fact valid for a broader class of metasurface response. The remaining parameters of dielectric loss and periodicity contribute significantly to the function of the metasurface. This work fixes the height, periodicity, and radius – thus the radiative rate – and varies the dielectric loss.

In the current understanding of the DMS a design pathway can be identified. First a working frequency and material for the resonator are identified. The choice of *ω*_{0} and material, i.e. $\tilde{\mathrm{\u03f5}}$, sets *r, h* by Eqs. (1) & (2). The remaining geometric parameter is periodicity (*a*), typically bounded by 2*r < a < λ*_{0} in order to maintain mode confinement and avoid significant diffraction effects, which sets the radiation rate of the structure. Finally, with the intended function of the metasurface in mind, a dielectric loss is chosen.

We perform S-parameter simulations of a unit cell of the dielectric array for silicon cylindrical resonators with different loss values and dimensions as noted in Fig. 1(a). Figure 1(c) and (d) shows the simulated transmission and absorption, respectively, as a function of frequency for a number of different dopings ranging from *n _{min}* = 1.275 × 10

^{13}to an optimal doping of

*n*= 7.374 × 10

_{opt}^{15}. Remarkably – as seen in Fig. 1(c) & (d) – a metasurface can be designed as a Huygens surface (blue curves) or, with the addition of a small amount of loss, the same metasurface can become an absorber (red curves). A maximum transmission value of 93.7% is obtained at 0.986 THz with an associated phase of 1.3

*π*radians for a doping value of

*n*, and a maximum absorption of 99.96% at 1.02 THz for

_{min}*n*.

_{opt}The experimentally realized metamaterial consists of an array of silicon cylinders arranged in a square lattice on a thin polydimethylsiloxane (PDMS) substrate (Fig. 1(a) & (b)). Four separate metamaterials were fabricated, with varied concentrations of dopant atoms to achieve distinct cases of dielectric loss – see Table 1. All DMS’s fabricated have the dimensions *a* = 172*µm*, *r* = 60*µm*, *h* = 50*µm*, on a 23*µm* PDMS substrate (Fig. 1(b)). The resistivity of the silicon substrate was measured for each of the four experimental cases (Table 1) allowing estimation of the dopant concentration [19], and subsequently calculation of the complex permittivity of each material. In the THz regime semiconductor permittivity is well described by a Drude model [20] as, $\tilde{\mathrm{\u03f5}}(\omega )={\mathrm{\u03f5}}_{\infty}-{\omega}_{p}^{2}/\omega (\omega +i{\omega}_{s})$ where, for silicon *ϵ*_{∞}= 11.7 is the frequency independent permittivity, *ω _{p}* is the plasma frequency, and

*ω*is the scattering frequency, defined as

_{s}*ω*= 1/

_{s}*τ*where

_{c}*τ*is the average collision time of electrons. The plasma frequency

_{c}*ω*is defined as,

_{p}*n*is the dopant concentration,

*e*is the fundamental charge,

*ϵ*

_{0}is the permittivity of free space, and

*m*

^{∗}is effective mass. The loss tangent was also calculated for the experimental cases and, incorporating the Drude permittivity, can be written as,

## 3. Experiment and simulation

With the origin of dielectric loss in the studied system presented in Eqs. (3) and (4), we experimentally explore its impact on the DMS system. The frequency dependent response of the DMS is simulated for the four cases given in Table 1, and fabricated samples are measured by THz time domain spectroscopy (TDS). The results (Fig. 2) show a varied scattering response, with the lowest loss achieving the lowest absorption (Fig. 2(a)), and moderate material loss achieving the highest absorption (Fig. 2(d)). Experimental results (Fig. 2 solid curves) agree well with simulated absorbing response (Fig. 2 dashed curves). The low loss (Fig. 2(a)) and high loss (Fig. 2(d)) results realize a shift between the transmission and reflection minima, and thus the metasurface is not strongly absorbing for these cases. Notably there is no significant shift in the peak absorbing frequency in the four cases.

An understanding of the effects of loss on the scattering properties of the DMS is furthered by simulation across the frequency and tan *δ* parameter space, and is supported by representation of the DMS in terms of its effective material constants of refractive index $(\tilde{n})$ and wave impedance $(\tilde{Z})$, which apply generally [21]. The scattering and material parameters are obtained from full wave electromagnetic simulations in CST Microwave Studio 2017, utilizing a frequency domain solver with unit cell boundary conditions. In Fig. 3 we present the index of refraction $(\tilde{n}={n}_{1}+i{n}_{2})$ and wave impedance $(\tilde{Z}={Z}_{1}+i{Z}_{2})$, which are defined at normal incidence for a system of thickness *d* = *h*. From simulation of the DMS’s scattering properties, the material parameters $(\tilde{n},\tilde{Z})$ can be obtained [22].

The extracted material parameters are presented as 2D colormaps in Fig. 3, with each vertical axis denoting the frequency in THz and horizontal axis the loss, i.e. tan *δ*. In Fig. 3(a) the real part of the wave impedance is shown, normalized by the free space value *Z*_{0} = 377Ω. The colorbar scale shown to the right of the plot details the values, and it can be observed that at a frequency of *ω*_{0} =1.02 THz (horizontal dashed line) and an optimal loss of tan *δ _{opt}* = 0.058 (vertical dashed line), we find a value of

*Z*

_{1}= 0.953

*Z*

_{0}– shown as the solid black contour curve in Fig. 3(a). In Fig. 3(b) the imaginary part of the refractive index is shown, normalized by its maximum value, and we find a value of ${n}_{2}=0.8577{n}_{2max}$ at

*ω*

_{0}and tan

*δ*. The absorption is plotted in Fig. 3(c) – a black contour line shows

_{opt}*A*= 99%, and we achieve a maximum value of 99.96%.

The description of the metamaterial in terms of its effective material constants further allows for an equivalent, simplified system to be considered. Two aspects of the system which offer insight are the material Q-factor, and field decay in a homogeneous material. In general the material Q-factor can be expressed as the inverse of a material’s loss tangent, with contributions from magnetic and electric response adding in parallel [23]. For the metasurface effective medium the total quality factor can be defined as,

where ${Q}_{E}^{-1}$ is the quality factor due to the electric response, and ${Q}_{M}^{-1}$ is the quality factor for the magnetic response. Both ${Q}_{E}^{-1}$ and ${Q}_{M}^{-1}$ can be expressed in terms of the effective optical constants as,Wave impedances shown in Eq. (6) are normalized by *Z*_{0}, and $\tilde{n}$ is the index of refraction. We thus find a total quality factor using Eq. (5) for our DMS of *Q _{T}* = 0.301 (Fig. 3(d)) at peak absorption. It is interesting to note that, for the simple harmonic oscillator, critical damping occurs at

*Q*= 0.5. In the considered case two modes exist which, if independent, would each be expected to attain critical damping. In a linear system Q-factors add inversely, leading to a total quality factor of 0.25 in the two oscillator system, which is comparable to the calculated value of

*Q*= 0.301.

_{T}We next consider the spatial dependence of the fields within an effective medium described by the index of refraction $\tilde{n}$ retrieved from S-parameter simulations. The fields of an electromagnetic wave propagating in the positive z direction may be written as,

*k*

_{0}=

*ω*

_{0}/

*c*,

*E*

_{0}is the amplitude of the electric field, and

*E*

_{0}/

*Z*the amplitude of the magnetic field in the homogeneous medium. We plot the spatial dependence of the fields for bulk silicon (Fig. 4(d)), and our effective material with the properties of the DMS (Fig. 4(a), (b)) for each value of tan

*δ*, as indicated in the legend. We plot field curves for four different loss cases, which we describe as: under-damping (tan

*δ*= 0.0001) and (tan

*δ*= 0.0253), over-damping (tan

*δ*= 0.25), and critical damping (tan

*δ*= 0.0582). In the critical damping state, perfect absorption is achieved with the normalized amplitude of the E-field decaying to a calculated value of 6 × 10

^{−3}, (from its value at

*z*= 0), at a distance of

*z*=

*h*, i.e. the height of the cylinder (Fig. 4(a), red curve).

The physics underlying the critical behavior of the system may be further explored using temporal coupled mode theory (TCMT) and eigenfrequency simulations (Fig. 5(b)). Figure 5(a) shows the absorptivity (A), reflectivity (R), and transmissivity (T), as a function of both doping (bottom axis) and tan *δ* (top axis). We find high *T* for the low loss case, and a maximum *A* at tan *δ* = 0.058. We also determined the loss dependence of the complex eigenfrequency $(\tilde{\omega}={\omega}_{1}+i{\omega}_{2})$ using the eigenvalue solver in COMSOL 5.3 for both the EH and HE modes. Eigenvalue simulations may be performed with or without material loss, thus yielding a dissipative damping rate (Δ) and radiative damping rate (Γ), where *ω*_{2} = Δ + Γ (Fig. 5(b)). [15] The eigenfrequency results complement a single frequency (1.02THz) analysis of the metasurface scattering (Fig. 5(a)), and three distinct regions are identified in Fig. 5 – under-damped (Δ/Γ < 1), absorbing (shaded), and over-damped (Δ/Γ > 1), with high absorption occurring at (Δ = Γ) for both the EH (red curve) and HE (blue curve) modes. Notably, both the EH and HE modes have a strong dependence on doping, with the quantity Δ/Γ increasing nearly linearly for both modes proportional to doping.

## 4. Discussion

The applicability of effective medium theory to the all-dielectric metasurface investigated is an important finding; the results of the effective material extraction and subsequent analysis are self-consistent, and offer powerful insight into the DMS system. Several factors elucidate the DMS response. We calculate a penetration depth from the refractive index for the critically damped case, $(\tilde{n}({\omega}_{0})=2.896+4.785i)$, and find a value of 9.782*µm* at perfect absorption. Further, the phase relationship between the electric and magnetic field may be found from the effective impedance ${\tilde{Z}}_{r}({\omega}_{0})=0.953+0.038i$, and remarkably a value of $arg(\tilde{Z})=0.011\pi $

rad = 1.98° is found. Thus *E* and *H* remain in-phase as they exponentially decay within the absorber – in stark contrast to the case for good conductors, where there is a *π/*4 phase shift, i.e. the magnetic field lags the electric field by 45°. These findings are consistent with the concept of an ideal absorber, and the in phase nature of the fields further verifies our result that the electric and magnetic dipoles are uncoupled, and dissipate energy equally. Indeed the total Q-factor obtained from the effective material parameters of the metasurface give a result suggesting that the DMSA acts as a set of uncoupled critically damped oscillators. Furthermore the excellent agreement between the calculations based on the effective homogeneous material (as seen in Fig. 4(a) & (b)), and the classical example of the damped harmonic oscillator strongly support the validity of this method. The effective material analysis also reveals a unique feature of the system, namely in the critically damped condition the electric and magnetic fields decay at identical rates (Fig. 4(a) & (b) – red curves). This finding offers a surprisingly simple description of a perfectly absorbing system, and its broader applicability to Huygens’ metasurfaces is of continued interest.

The dissipative rate (Δ) of the DMS is proportional to the doping density *n*, and determined by Eqs. (3) and (4). The even mode, also termed the electric mode, clearly owes its loss to the complex permittivity; however, the odd magnetic mode is less obviously related to *ω _{s}*. Careful examination of the fields [10, 18] reveals a characteristic circulating electric field, therefore the loss due to the magnetic mode also results from the complex permittivity, and thus depends on

*n*. Thus the mechanism of high absorption in the DMS is due to the balance between dissipative and radiative rates, and is tuned directly by doping as illustrated in Fig. 5(b), where the ratio of Γ to Δ is proportional to doping concentration over most of the parameter space. This finding results in a direct design pathway as the data suggest that Δ ∝

*n*, with the function of the DMS in the low loss regime varying smoothly with Δ.

The versatility of the all-dielectric system is summarized well in Fig. 1(c) & (d) and Fig. 5(a), where the smooth transition from a Huygens surface to an absorber is evident. In the low loss regime, where high (≈ 78%) transmission occurs with low (≈ 0%) levels of absorption (Fig. 5(a)), the system is dominated by the radiative rate (Fig. 5(b), Δ/Γ < 1). Importantly in this low loss regime, significant phase modulation can be achieved – a key requirement for the Huygens surface. High absorption is achievable (Fig. 1(d), Fig. 2(b), and Fig. 5(a) – gray region); however, absorption does not increase monotonically with material loss. Instead, absorption increases rapidly from low loss, and trails off slowly after peak value. This is the expected result, as the balance of the radiative and dissipative effects for a fixed geometry occurs at a distinct value of material loss [15, 17], and after achieving perfect absorption, the system is dominated by dissipation (Fig. 5(b)). Large phase variation, coupled with the dependence on loss tangent, indicate that phase and function can be tuned by adjusting carrier concentration, highlighting the utility of semiconductors as a base material for dynamic, multi-functional metasurfaces.

The role of loss in the system can be summarized in a straightforward manner – loss determines function. The radiative rate of the system is fixed by the geometry given an operational wavelength. This radiative rate, once set, is counterbalanced by a dissipative rate, which is the sum of all non-scattering loss processes. In the all-dielectric system increasing dielectric loss directly increases the dissipative rate, thus the low loss regime is a radiative system, where the amplitude and phase of the output of the system can be set by choosing some small value of loss. As loss is increased the system becomes balanced, and achieves a critical state when Γ = Δ, where perfect absorption is achieved. The remainder of the parameter space is primarily dissipative in nature, with a slow relaxation of the matching condition.

## 5. Conclusion

The present work has demonstrated that the function of all-dielectric metasurfaces can be controlled by altering the dielectric loss of the base materials. Semiconductors are an ideal material to fashion a dynamic all-dielectric metasurface and offer tremendous versatility, as the system can be smoothly transitioned from functioning as a low-loss Huygens surface, to a perfect absorber. This diverse set of possible responses has been computationally explored and experimentally verified. Furthermore an effective material description is developed and illustrates the system’s behavior as a harmonic oscillator, with critical damping at perfect absorption. This analogy, and accompanying analysis of the Q-factor of effective materials offers powerful insight into absorbing systems, and provides an interesting avenue for future study. Finally, the extreme versatility exhibited by the metasurface is a unique feature of this class of all-dielectric metamaterial, and offers exciting possibilities for novel, multi-functional and dynamic metamaterial systems.

## 6. Experimental design

#### 6.1. Fabrication

The absorbers were fabricated by identical processes on four silicon-on-insulator (SOI) wafers of varied boron dopant concentration. The device layer was patterned using a Bosch process to achieve the necessary etching depth and straight sidewalls. The cylinders were released from the SOI wafer by etching the buried oxide with 49% hydrofluoric acid to achieve a suitable undercut. The device is then flipped and bonded to PDMS on glass substrate by curing at 100°*C* for 1 hour. Subsequent application of a lateral force results in the release of the cylinders from the SOI wafer, leaving an array of cylinders on a thin (23*µm*) PDMS substrate. Finally the device is bonded to a free standing, thick (3*mm*) PDMS support substrate, and released from the glass by dissolving a sacrificial layer of photo-resist in acetone.

#### 6.2. Material characterization

The DMS’s were characterized via THz time-domain spectroscopy. The devices were placed in the focal plane of the optical path in both transmission and reflection. Diffuse scattering is assumed to be negligible, giving *A* = 1 −*T* − *R*. The measurements were taken in a dry air purged environment. All reported results are normalized measurements with a gold mirror as a reference for reflection, and an open channel (empty sample holder in a dry air purged environment) as a reference for transmission. In order to minimize the effects of detector alignment between measurements each reported curve is the averaged result of three measurements alternating between reflection and transmission. The resistivity of the silicon was measured on device areas prior to patterning using standard 4-point probe techniques.

In the THz measurements, a deviation in the absorption spectrum occurs near 1.1 THz in some spectra, and is centered between two atmospheric absorption lines. These atmospheric effects likely contribute; however, variation in the surface normal with respect to the incident THz beam is expected to produce a small angle dependent peak, spatially summing to the feature observed. We attribute differences in simulated and experimental absorption to variations in the surface normal of our samples, due to the PDMS substrate which lacks rigidity. Additionally any off-normal scattering processes enter into the absorption spectrum as these effects would not be captured at the detector.

## Funding

Department of Energy (DOE) DE-SC0014372. 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. **D. R. Smith, W. J. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, “Composite Medium with Simultaneously Negative Permeability and Permittivity,” Phys. Rev. Lett. **84**, 4184–4187 (2000). [CrossRef] [PubMed]

**2. **D. Schurig, J. J. Mock, B. J. Justice, S. a. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science **314**, 977–980 (2006). [CrossRef] [PubMed]

**3. **N. I. Landy, S. Sajuyigbe, J. J. Mock, D. R. Smith, and W. J. Padilla, “Perfect metamaterial absorber,” Phys. Rev. Lett. **100**, 1–4 (2008). [CrossRef]

**4. **T. Maier and H. Brückl, “Wavelength-tunable microbolometers with metamaterial absorbers,” Opt. Lett. **34**, 3012–3014 (2009). [CrossRef] [PubMed]

**5. **N. Liu, M. Mesch, T. Weiss, M. Hentschel, and H. Giessen, “Infrared perfect absorber and its application as plasmonic sensor,” Nano Lett. **10**, 2342–2348 (2010). [CrossRef] [PubMed]

**6. **C. M. Watts, X. Liu, and W. J. Padilla, “Metamaterial electromagnetic wave absorbers,” Adv. Mat. **24**23 (2012).

**7. **C. Shemelya, D. Demeo, N. P. Latham, X. Wu, C. Bingham, W. Padilla, and T. E. Vandervelde, “Stable high temperature metamaterial emitters for thermophotovoltaic applications,” Appl. Phys. Lett. **104**201113 (2014). [CrossRef]

**8. **J. A. Schuller, R. Zia, T. Taubner, and M. L. Brongersma, “Dielectric metamaterials based on electric and magnetic resonances of silicon carbide particles,” Phys. Rev. Lett. **99**, 1–4 (2007). [CrossRef]

**9. **S. Jahani and Z. Jacob, “All-dielectric metamaterials,” Nature **11**, 23–36 (2016).

**10. **X. Liu, K. Fan, I. V. Shadrivov, and W. J. Padilla, “Experimental realization of a terahertz all-dielectric metasurface absorber,” Op. Exp. **25**, 191 (2017). [CrossRef]

**11. **I. Staude, M. Decker, E. Rusak, D. N. Neshev, and I. Brener, “Active Tuning of All-Dielectric Metasurfaces,” ACS Nano **9**, 4308–4315 (2015). [CrossRef] [PubMed]

**12. **P. Gutruf, C. Zou, W. Withayachumnankul, M. Bhaskaran, S. Sriram, and C. Fumeaux, “Mechanically tunable dielectric resonator metasurfaces at visible frequencies,” ACS Nano **10**, 133–141 (2016). [CrossRef]

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

**14. **E. Kallos, I. Chremmos, and V. Yannopapas, “Resonance Properties of Optical All-Dielectric Metamaterials Using Two-Dimensional Multipole Expansion,” Phys. Rev. B **86**, 245108 (2012). [CrossRef]

**15. **X. Ming, X. Liu, L. Sun, and W. J. Padilla, “Degenerate critical coupling in all-dielectric metasurface absorbers,” Opt. Exp. **25**, 24658 (2017). [CrossRef]

**16. **B. Neuner, D. Korobkin, C. Fietz, D. Carole, G. Ferro, and G. Shvets, “Critically coupled surface phonon-polariton excitation in silicon carbide,” Opt. Lett. **34**, 2667–2669 (2009). [CrossRef] [PubMed]

**17. **C. Qu, S. Ma, J. Hao, M. Qiu, X. Li, S. Xiao, Z. Miao, N. Dai, Q. He, S. Sun, Y. Zhang, and L. Zhou, “Tailor the functionalities of metasurfaces based on a complete phase diagram,” Proceedings of IEEE International Workshop on Electromagnetics , **235503**, 1–6 (2016).

**18. **K. Fan, J. Y. Suen, X. Liu, and W. J. Padilla, “All-dielectric metasurface absorbers for uncooled terahertz imaging,” Optica **4**, 601 (2017). [CrossRef]

**19. **B. L. Anderson and R. L. Anderson, *Fundamentals of Semiconductor Devices* (McGraw-Hill College, 2004), 1st ed.

**20. **M. van Exter and D. Grischkowsky, “Optical and Elelctronic Properties od Doped Silicon from 0.1 to 2 Thz,” Appl. Phys. Lett. **56**, 1694–1696 (1990). [CrossRef]

**21. **C. R. Simovski, “Material parameters of metamaterials (a Review),” Opt. and Spect. **107**, 726–753 (2009). [CrossRef]

**22. **D. R. Smith, D. C. Vier, and C. M. Soukoulis, “Electromagnetic parameter retrieval from inhomogeneous metamaterials,” Phy. Rev. E **71**, 036617 (2005). [CrossRef]

**23. **R. Simon, J. Whinnery, and T. Van Duzer, *Fields and Waves in Communication Electronics* (John Wiley & Sons, 1994), 3rd ed.