A theoretical investigation on reciprocity-inspired wide-angle spectrally-selective THz absorbers augmented by anisotropic metamaterials

In this paper, a theoretical framework relying on the reciprocity theorem is proposed to accurately design a spectrally-selective THz superstrate-loaded metamaterial absorber (SLMA) exhibiting wide-angle feature. By leveraging high-order Floquet harmonics in a generalized transmission line model characterizing the conventional metamaterial absorbers (MAs), it is demonstrated that MAs suffer from impedance mismatch, especially at near grazing angles. From an impedance matching viewpoint, this major challenge is tackled in this paper via two different designs, exploiting a magneto-electric anisotropic Huygens' metamaterial and a multilayer dielectric structure at a certain distance over the MA plane. The numerical results corroborate well the theoretical predictions, elucidating that the proposed SLMA significantly broadens the angular performance of the MA up to near grazing angles (about 80°), where high absorptivity is still achieved in both principal planes. The deteriorating effect of diffraction modes has been comprehensively analyzed. In comparison to the previous wide-angle MA reports based on intricate particle geometries and brute-force optimizations, the proposed design features a straightforward semi-analytical algorithm, which can also be re-developed for microwave, mid-infrared, and optical frequency bands and for any type of MA element. The proposed SLMA would be very promising for various wavelength-selective applications such as sensors and imaging.

www.nature.com/scientificreports/ property of the MA is solely dependent on the rotational symmetry of the constituent units 23,24 . Acquiring perfect absorption with the angle-insensitive privilege, however, is a much more daunting task 19 due to the wave-MA impedance mismatch appearing near the grazing angle where the reflectance quickly approaches 100%. Thus, at extreme angles, matching with conventional methods becomes difficult if not impossible, so that for existing wide-angle MAs, the region close to the end-fire direction is difficult to cover. Recently, an enormous amount of effort has been made to enhance the angular features of MAs with engineered geometries, e.g. split-ring-cross resonator 25 , four-fold rotational symmetric electric resonator with a cross-printed bottom 26 , deep subwavelength unit cell in a multilayer 27 , surrounding via array 28 , and circular sector 29 . Although the sensitivity to the polarization and angle of incidences is reduced to some extent, the vast majority of reports in all spectra deal with sophisticated, optimization-based, and costly designs, as they require the deformation of all array elements [30][31][32][33][34][35][36][37] .
Further, these approaches only support certain types of array elements, probably creating features which make them difficult for applications over a large area 38 . For instance, based on an anisotropic perfectly impedancematched negative-index material, a metamaterial-based approach was described 39 for making a wide-angle absorber of infrared radiation. Nevertheless, due to the single layer nature of the absorber, it operated only for s-polarized (single-polarization) oblique incidences up to 45°. Besides, via the presented approach, a lossy anisotropic metamaterial must be designed to expose specified real and imaginary parts of the permittivity and permeability at the operating frequency and, thus, the metamaterial design would be very complex. Shen et al. 40 reported a compact single unit cell consisting of three nested electric closed-ring resonators as a microwave tripleband absorber. However, as reported by them, the designed absorber was valid to a wide range of incident angles up to only 50° for both TE and TM polarizations. Meanwhile, the design geometry was selected heuristically, and no exact analysis was presented to formulate (or at least justify) the wide-angle performance of the designed MA. Finally, it worked at the microwave frequency band, not THz spectrum. A simple metamaterial-based wide-angle plasmonic absorber was introduced, fabricated, and experimentally characterized 41 using angle-resolved infrared spectroscopy. The main focus of that paper was not on expanding the angular operating range of the absorber, but on low-cost designing of a plasmonic absorber via nano-imprint lithography for manufacturing large-area samples and so, the angular performance of the absorber was not high enough for diverse applications. More importantly, the design parameters were obtained based on numerical optimizations, and the semi-analytical formulas given in that manuscript were only provided for several justifications about the impedance matching.
There are also other studies that have addressed wide-angle metamaterial absorbers at different parts of the spectrum [42][43][44] . Nevertheless, in most of them, the design of absorbing metamaterial unit cells relied on trial and error steps, no rigorous theoretical formulation supported the study, and the design was based on a heuristic approach. Another possible solution to realize artificial magnetic response is to use a high-impedance surface (like mushroom structure) 45,46 . Most of these designs are for the microwave frequency range, where each patch is connected to the ground plane with a thin conducting via, forming a "mushroom". The presence of vias is important for the operation at oblique incidence, where they become useful also in absorber applications. Olli Luukkonen et al. 45 designed a mushroom-like absorber with a high-permittivity substrate (obviously with highcost fabrication) whose absorption efficiency remained high for incidence wave angles up to 60° at microwave frequencies.
In the case of a low-permittivity substrate (logical manufacturing cost), the designed absorber failed to achieve satisfactory results for TE polarization at 60°4 5 . Consequently, there is a critical need to develop a new architecture revealing wide-angle THz absorptivity combined with simplicity and feasibility. However, little attention has been paid to discussions about the theories or mechanisms of wide-angle absorption, which is very crucial for potential research on MAs.
Here, a combined version of the wide-angle matching layer in phased-array antennas 47 , and reciprocity theorem 48 from classical electromagnetics is theoretically established to address a spectrally selective semiomnidirectional THz absorber. We will serve magneto-electric anisotropy as an additional degree of freedom to provide the required wave-absorber impedance matching at near grazing angles. In a more simple design with respect to the perfect match layers, we require a loss-less anisotropic metamaterial exhibiting only the predetermined real parts of the permittivity and permeability at the operating frequency. Indeed, we present two matching superstrate layers with uniaxial tensors located at a certain distance from a pre-designed MA, a simple configuration physically accomplishable via Huygens' metamaterial or multilayer dielectrics. For demonstration purposes, the designed wide-angle matching layer is located at a certain distance from a printed-dipole MA to provide minimum reflection levels for a wide angular range of oblique illuminations. We will demonstrate that, with the appropriate selection of the transversal and longitudinal anisotropic parameters characterizing the superstrate layer, strong absorptivity can be acquired at 2 THz upon illumination by a wide range of incident angles (up to 80°). The superstrate synthesis process is aided by a parameter extraction technique, accurately extracting the effective uniaxial components. The proposed superstrate-loaded MA (SLMA) does not require a complex design, and more importantly, can be applied to any type of MA and sensor element.

Results
Analytical scheme. In this section, without loss of generality, we will exemplify the proposed design method by leveraging an infinite array of spectrally-selective dipoles as the basic building block of the MA. With the expansion of printed technologies, dipole arrays are known as cost-effective, lightweight, and low-profile meta-devices 49 , so that the advantages of planar MAs becomes much more tangible when employing dipole elements. Furthermore, the use of dipole elements allows an accurate analytical solution for predicting the THz response of MA augmented by auxiliary anisotropic layers. Figure 1 indicates the periodic MA array comprising of zero-thickness dipoles with the length of l, width of w , and periodicity of p x and p y along x-and y-directions, respectively. The array elements are located on a polystyrene foam 50 (n d = 1) with the thickness of h 1 terminated by a gold ground plane ( σ = 4.56 × 10 7 s/m ) to block the transmission. A finite-size ultra-thin sheet of ran- www.nature.com/scientificreports/ dom silver nanowire (AgNw) network 51 responsible for dissipating the captured energy terminates the middle of each dipole element. The AgNw layer can be thought of as a resistive load for the array of dipoles, which is characterized by its surface resistivity, R s (Ω/m 2 ). In a strict sense, the equivalent lumped resistance of the lossy AgNw depends on the unit cell shape such that the scattering surface area per unit is larger in the uniform sheet model than the physical scatterer. As an acceptable estimate 52 , the addition of lossy AgNw layer can be circuitally modelled through a termination resistance of R L = R s p 2 x /S wherein S refers to the surface of the lossy layer along the direction of the induced current. Thanks to the reciprocity theorem, the array of dipoles with the best transmission specifications in the radiating mode will similarly exhibit the best absorption performance when it is suitably terminated in the receiving mode 53 . Consequently, to determine the angular-dependency of the reflectivity for dipole elements in the absorbing mode, one can utilize the scan impedance of the dipoles in the transmitting mode which is defined as the observed impedance when looking into the terminals of a driven dipole element surrounded by an infinite array of similarly driven elements 54 . In this case, the periodic modality of the array allows us to describe the near-field interactions of the dipole array via the Fourier series expansion of the radiating fields. In this way, the TM-polarized fields can be expressed as the superposition of the fields radiated by the y-directed dipole surface currents and the image fields caused by the ground plane in plane waves indexed by different harmonics in which, η and k refer to the characteristic impedance and the wavenumber of the host medium, respectively; and k xm = k sin θ cos ϕ + 2mπ/p x , k yn = k sin θ sin ϕ + 2nπ/p y ; and The surface current density of each unit cell is zero except on the dipole element, which, in turn, can be described by a complex Fourier series. Assuming a sinusoidal loop of current with the maximum (I 0 ) at the dipole center, K mn coefficients can be subsequently disclosed using the orthogonality principle 4 On the other hand, the scan impedance of each radiating dipole is calculated by sin k x,m w/2 k x,m w/2 cos k y,n l/2 www.nature.com/scientificreports/ which yields The above sum can be immediately interpreted as the contribution of the propagating/non-propagating Floquet harmonics indexed by m and n to the input impedance of the dipole array. Here, the last term represents the additional contribution of the reflections off the ground plane and 55 where cosθ x = k x /k and cosθ y = k y /k. In the receiving/absorbing mode, the MA must be suitably terminated by a resistive load (A piece of AgNw lossy sheet is utilized in this paper) to harvest the maximum power impinging on the array. It should be noted that another lossy material in the THz spectrum, such as graphene 9 , could also be utilized to dissipate the trapped power. The reciprocity condition ensures that the reflectivity seen by a plane wave looking into the MA array (in the absorbing mode) is equivalent to the reflectivity seen looking into the terminal of the dipoles (in the transmit mode) in which R L denotes the terminated/source impedance. The absorption spectrum of the dipole MA can be calculated from A = 1 − R − T, in which R =|S 11 | 2 and T =|S 12 | 2 are the reflectance and transmittance of the array, respectively. Nevertheless, since the gold ground plane thickness is much larger than THz skin-depth, we have no transmission (T = 0) across the frequency band of interest. The related parameters of the designed MA are p x = p y = 0.5 , l = 0.48 , w = 0.02 , h 1 = 0.25 , and R L = 50 Ω where denotes the operating wavelength corresponding to f = 2 THz. A MATLAB program was prepared to implement the input (scan) impedance of the MA array whose accuracy with the involved multiple harmonics is compared with the benchmark full-wave numerical solution (CST commercial software). An adequate number of Floquet harmonics, at least including m, n ∈ {0, ±1, ±2, ±3, . . . , ±10} , is considered here to reach acceptable accuracy in the analysis, as the other harmonics were found to have negligible contributions. On the other hand, in the full-wave modelling of the proposed MA, periodic boundary conditions are applied in both the x-and y-directions while the Floquet port illuminates the structure along the z-direction. Figures 2a, b compare, respectively, the analytical and numerical results of absorptivity for the employed MA in both E-(yOz plane) and H-planes (xOz plane). As the inset of these figures show, an acceptable agreement is noticed between the reference solutions and the analytical results. The E-plane absorptivity deteriorates at the boresight direction since the scan impedance of the array is high in this angular regime (about 150 Ω). Besides, even though the absorptivity variation of the array is slightly smoother in the E-plane, the angular dispersion is not weak in both planes. This observation can be attributed to the mutual coupling between the elements where the reactive component of the input impedance for the employed MA increases with the scan angle 49 . Therefore, due to the air-absorber impedance mismatch, the MA array of Fig. 2 exhibits an angle-sensitive performance in both principal planes, where reflectivity approaches 100% near the grazing angle.

The near-field interactions between MA array and uniaxial superstrate layer
For optimal MA performance upon illumination by oblique incidences, it is highly desirable to counteract the strong angular dependency of the input impedance of the dipole MA with an auxiliary anisotropic layer located at a certain distance, h 2 , above the array (see Fig. 3a). Assuming an anisotropic superstrate layer characterized by diagonal (uniaxial) magnetic and dielectric tensors 56 www.nature.com/scientificreports/ will substantially simplify the synthesis and enhance the feasibility of our design. Since the evanescent Floquet harmonics in the vicinity of the array and the superstrate layer can highly interact with the discontinuity caused by the anisotropic layer, the resultant effect on the scan impedance of the MA becomes greater as the superstrate approaches the array. Inspired by the form of Eq. (6), a network of uncoupled transmission lines can be served to take into account such near-field interactions between the MA elements and the anisotropic superstrate layer (with the thickness t above the array). Keeping the reciprocity theorem in mind, the fields of the dipole array (in the radiation mode) can be expanded with a set of distinct Floquet harmonics, each of which can be represented by a simple equivalent transmission line model. It is for the first time that the role of evanescent Floquet harmonics is theoretically assessed when designing MAs augmented by auxiliary anisotropic layers. Here, considering different media above the dipole elements, a certain number of clustered cascaded lines with different characteristics describe each harmonic (Fig. 3b). According to the developed semi-analytical relations for MA (array of dipoles terminated by the lossy AgNw sheet), it can be electromagnetically incorporated before the uncoupled transmission lines modelling the overlaying media at each Floquet harmonic. Upon being affected by the anisotropic impedance-matching layer, the new scan impedance at the input aperture discontinuity of the array can be calculated as a series summation of the individual scan impedances belonging to each Floquet harmonic (Fig. 3b). The latter case can be simply computed based on several transmission line impedance transformations  www.nature.com/scientificreports/ to move the impedance at the boundary z = h 1 + h 2 + t to that of z = h 1 . Here, we will treat the problem in separate cases of TM and TE polarizations where only the ε y , ε z , and µ x components matter for TM polarization, while only the ε x , µ y , and µ z components play the matching role for TE polarization. In this fashion, the EM properties of the anisotropic layer, i.e. the wave impedance and normal component of the wave vector, can be written as 47 Moreover, the addition of the anisotropic layer also re-defines Eq. (6) considering the multiple reflections occurring between the nearby interfaces 49 in which The impedances for each Floquet harmonic, η TE mn and η TM mn , are given by k/ν mn and ν mn /k , respectively, and Ŵ TE/TM out,mn stands for the reflection observed looking into the anisotropic layer. Besides, Z TE/TM b represents the transferred impedance of the free space medium at the boundary ( z = h 1 + h 2 + t) to that of z = h 1 + h 2 . Hereafter, Eq. (12) can serve to compute the modified scan impedance and the absorptivity of the proposed SLMA, i.e. the new absorber architecture augmented by the uniaxial anisotropic superstrate. As the dimensions of the proposed meta-atom are comparable with the wavelength (p = 0.5λ), care should be taken about the presence of high-order Floquet harmonics (diffraction modes) in the calculation of absorptivity values, especially in the oblique incidence situation. A comprehensive discussion is given in Supplementary Information A. In this case, the absorptivity of the proposed structure can be computed based on: in which, the second term denotes the co-pol, and the third one refers to the cross-coupling impacts of the high-order Floquet harmonics. Also, parameters m and n denote the index of each diffraction mode. All results presented throughout this paper have been computed considering a truncated number of diffraction modes.
In the first step, we prepared a constrained optimization to seek optimum parameters revealing the best wideangle absorption performances among all the possible solutions that may, however, be arduous to be realized. The magneto-electric parameters of the anisotropic matching layer are allowed to be freely varied within intervals, excluding extreme values. The constraints used in the optimization at f = 2 THz are tabulated in Table 1. This optimum solution is utilized only as a benchmark to demonstrate the wide-angle impedance matching caused by the anisotropic superstrate.
For a set of 300 randomly generated starting points, the local minima of the cost function are repeatedly sought. Figure 4a, b plot the resultant absorption curves versus the angle of incidence in both E-and H-planes, respectively. Evidently, the theoretical and numerical results perfectly overlap each other, confirming the validity of our proposal for analyzing the proposed SLMA. Consequently, the analysis yields exact results,

Realizing the wide-angle SLMA
The anisotropic parameters used for the results of Fig. 4a, b are best-case solutions assuming that all dielectrics with arbitrary EM parameters are available, which is not true in real-life scenarios. To circumvent this unrealistic assumption and due to the importance of practical considerations, this study discloses two distinct accomplishable proposals, i.e. electric-only and magneto-electric uniaxial superstrate layers, to implement the designed SLMA. Figure 5a, b illustrate the topology of two configurations of SLMA whereby the angular performance of the superstrate layers is computationally verified. To realize the electric-only uniaxial layer, we adopt the superstrate by periodically stacking two commercially available homogeneous materials with different permitivities and thicknesses of (ε r1 , ε r2 ) and (d 1 , d 2 ) along the perpendicular direction (see Fig. 5a). The equivalent tensor parameters of the multilayer dielectric system are retrieved in the absence of the array of dipole resonators. This is a common practice in analyzing the anisotropic superstrates loading a periodic array 47 . Based on the effective medium theory 57,58 , the equivalent tensor components of the dielectric permittivity can be derived as   www.nature.com/scientificreports/ On the other hand, to realize the magneto-electric anisotropic layer, we employ two parallel gold rings embedded in a host dielectric medium (polyethylene-high density with the refractive index of n = 1.54 59 ), separated by an air hole. The Huygens' metamaterial will strongly interact with both electric and magnetic incident fields, and its working principle can be individually evaluated for the tangential and normal components. The magnetic and electric parameters can be controlled by tuning the geometric parameters in a coupled manner 60 . In this case, the simple volumetric averaging of Eqs. (16a, b) is no longer applicable, and a more complicated approaches should be employed. Kramers-Kronig limitation is applicable to any material parameter intended to have a physical sense of a generalized susceptibility. Care should be taken that the selected system of describing parameters may fail to obey the relevant restrictions when its parameters, such as frequency-dispersive permittivity and permeability, are not sufficient to describe the actual response of the structure. In this case, a more extended theoretical description is required, taking spatial dispersion into account [61][62][63][64][65][66][67] .

numerical simulations
In this section, we aim to investigate the broad angular performance (or weak angular dispersion) of the SLMA realized by the multilayer dielectrics and the Huygens' metamaterial over the aperture plane of the MA. It should be noted that, although the Huygens' metamaterial is established to realize the best parameters given in Table 1, these parameters cannot be implemented by the multilayer dielectric system as it does not provide any magnetic response. Therefore, a new optimization procedure was employed, and the optimum dielectric tensor parameters were accordingly achieved. The best parameters were searched within a feasible database of THz materials available in practice, i.e., among the anisotropic parameters which are neither extremely high nor near-zero nor negative, as strong frequency dispersion and non-negligible losses will appear in these regimes. Finally, a specifically designed stack of Al 2 O 3 and PMMA layers was chosen to realize the required anisotropy of the superstrate. The optimum parameters are given in Table 2. Moreover, to find the optimal geometries of the Huygens' metamaterial mimicking the parameters of Table 1, extensive parametric sweeps were performed at the operating frequency of 2 THz. Being located far from the structure resonance, the proposed meta-atom almost treats in a non-dispersive routine in frequency. The main reason can be attributed to the fact that the magnetic polarization of a subwavelength metal ring has a flat frequency trend 47 . This can be mainly attributed to the fact that the metamaterial unit has no vias and operates far from any structural resonance. For further validation, the full-wave simulation results, as well as the theoretical predictions for the best solution of both SLMA designs, www.nature.com/scientificreports/ are displayed in Fig. 6a-d. The configuration of the simulation setup is shown in Fig. 5c in which the realized superstrate (double gold ring structure) and the dipole array are considered in a single simulation system, and the absorptivity of the whole system is extracted. As can be seen, nine sub-wavelength MM elements are located in each cell of the designed SLMA. The comparison depicts an almost perfect agreement between full-wave and semi-analytical results on the E-and H-planes, respectively. The results demonstrate that the magneto-electric Huygens' layer modifies the angular performances of the dipole MA array almost in a similar fashion to the multilayer superstrate for the E-plane, whereas due to providing both electric and magnetic behaviors, better response in the absorption has been observed on the H-plane. In this plane, the electric-only multilayer superstrate functions relatively well at lower scan angles, but its performance degrades in middle and higher scan angles. Overall, via the proposed designs, the absorbed power of the designed SLMA structures is dramatically increased in both  www.nature.com/scientificreports/ E-and H-planes when excited by a broad range of oblique incidences, a great capability that has not been analytically reported yet. Approximately, a 70% portion of the absorbed power is observed in the E-plane for oblique incidences up to 80° at the operating frequency. It is clear that our design is normalized to the wavelength and hence applicable to any range of the EM spectrum. However, one should be very careful with the model used for the metallization at THz as the nominal conductivity does not fit properly. Therefore, the practical aspects of the design can be further considered by reducing the nominal DC conductivity by order of magnitude to account for extra loss introduced by roughness. The corresponding results can be found in Supplementary Information B. To evaluate the sensitivity of the SLMA absorption (realized by Huygens' metamaterial) upon deviations in the anisotropic parameters of the employed superstrate (possible fabrication tolerances), we present the results in Fig. 7a, b for the performance of the optimized SLMA of Fig. 5b when the magnetic or/and electric tensor parameters deviate by 10% from the optimized solutions. Indeed, a 10% variation is considered for all the electric and magnetic tensor parameters. As the inset of this figure illustrates, the design is almost robust against the possible fabrication tolerances which may occur during the manufacturing processes, especially in the E-plane. Besides, since the frequency bandwidth is also a critical figure of merit in most THz systems, it is important to consider the implications of the anisotropic layers, whose particles' individual response may experience a rapid spectral change on the bandwidth of the dipole MA array. As discussed before, for the employed subwavelength particle realizing the Huygens' anisotropic metamaterial, the resonance frequency is far enough from the working frequency of 2 THz, so that the SLMA operates in a less frequency dispersive and less lossy regime of the anisotropic layer. In order to clearly reveal the spectrally-selective and wide-angle absorption properties of the proposed SLMA (realized by Huygens' metamaterial), we illustrate the absorption spectra at different angles of incidence in Fig. 8a, b. On the one hand, the value of the absorption peak is 0.98 for normal incidence. With increasing the incident wave angle, the value of the absorption peak remains interestingly high (around unitary) in both E-and H-planes. However, there is a slight frequency shift with an increase of θ, thus exhibiting minimal dispersion within the operating bandwidth. Obviously, the frequency bandwidth of the SLMA, which is almost dictated by the dipole elements, is not deemed as the main topic in this paper; nevertheless, the other engineered shapes for the primary MA can be selected to broaden the absorption band or for dual-polarized operation. In these cases, the EM field excited by the basic array can still be taken into account with semi-analytical models, when regular geometries such as rectangular and circular apertures and printed slots are involved. In addition, resorting to full-wave simulations is an alternative approach, which allows the analysis of arbitrarily shaped MA units 68 . Indeed, the proposals given in this paper do not necessarily postulate the primary knowledge about the analytical investigation of the MA element, as they can be numerically applied to any type of MA geometry.

Discussion
A wide-angle and spectrally-selective THz SLMA based on a fully analytical framework involving the near-field interactions between the Floquet harmonics and the superstrate layer has been presented in this paper for the first time. Using transmission line networks and relying on the reciprocity theorem, two designs based on an array of dipoles armed with uniaxial superstrate were sketched and analyzed. The multilayer dielectric system and the Huygens' metamaterial are separately utilized for realizing the electric-only and magneto-electric anisotropic superstrates, respectively. The comparison of the full-wave simulations and theoretical results for both designs revealed perfect similarity, verifying the proposed analytical scheme. It was demonstrated that, when suitably selecting the anisotropic parameters, the efficiency of MA could be substantially improved for a wide range of incident angles, including near grazing angles. As critical figures of merit, the effects of the frequency variation and the possible fabrication tolerances were investigated for the proposed SLMA. The proposals discussed in this www.nature.com/scientificreports/ paper do not necessarily postulate the primary knowledge about the analytical investigation of the MA element, as they can be numerically applied to any type of MA geometry when other features like broader bandwidths are highly demanded. Overall, the novelty of this paper was not only the design of a wide-angle THz MA but the proposal of a fully theoretical algorithm to expand the angular operation of a pre-designed absorber, regardless of its operating frequency. The main contributions of this paper are: (1) for the first time, our study accurately explored the interaction between the superstrate structure and the dipole array by incorporating the higher-order Floquet modes in the transmission line model; (2) the design is independent of the center frequency and can be re-performed for all spectra from microwave to visible; (3) the proposals given in this paper do not necessarily postulate the primary knowledge about the analytical investigation of the MA element, as they can be numerically applied to any types of MA geometry, e.g. to broaden the absorption band. The wide-angle property of the proposed SLMA would be very promising as an absorbing element for various applications such as THz imaging. Other appealing features can also be envisioned for the proposed structure, such as tunable functionality provided by the graphene 69-72 material.