This work focuses on finite element modeling (FEM) of a three-dimensional metamaterial used as an absorber for terahertz energy harvesting. The metamaterial consists of patterned pillars of an SU-8 dielectric photoresist coupled to a copper metal overlayer. Our study shows that the electromagnetic performance of the metamaterial is dependent on the following characteristic design parameters of the SU-8 dielectric: pillar height, bottom side length, and spacing between adjacent pillars. Using FEM, the metamaterial geometry is successfully optimized and the surface plasmon can be tuned to a peak frequency of 1.2 THz and a maximum terahertz absorption amplitude of 30%.

## I. INTRODUCTION

Metamaterials (MMs) are a type of fabricated synthetic material that process the ability to manipulate electromagnetic waves using its subwavelength repeating structures.^{1–4} Based on metamaterials’ structure, there are typically 1D, 2D, and 3D metamaterials. 1D metamaterials typically consist of planar layers of metals and dielectrics, which are commonly used in optics in the forms of filters or Bragg reflectors^{5} with an effective permittivity approaching zero. 2D metamaterials, which consist of periodic patterns of metal and dielectric, can be used to create materials with a negative index of refraction.^{6,7} The feature size of 1D and 2D metamaterials ranges from a few nanometers to a few hundreds of nanometers. If the targeted wave is on the scale of micrometers to millimeters, 3D structures are required to capture and localize it.^{8–10}

In terms of energy harvesting, the direct output from energy conversion is always an alternating voltage, so a metal-insulator-metal (MIM) diode^{11–13} is used in this study to convert alternating current (AC) to direct current (DC) for electrical loads. Besides, an energy harvesting device has to be closely tuned to its power source in order to achieve high conversion efficiencies.^{14} Metamaterials present an innovative conversion media that can selectively absorb the desired electromagnetic radiation with coherent waves.

The goal of this study is to use metamaterials as an absorber to harvest industrial waste heat. When an electromagnetic wave, in the form of heat or light, irradiates a metamaterial with a resonant frequency, an intense “plasma” of lower-frequency surface plasmons is generated at the metal/dielectric interface. Since surface plasmons are a type of free electron oscillation, their waves still exhibit electromagnetic characteristics. By optimizing the geometries of the metamaterial, the wide radiation spectrum can be selectively absorbed into a narrow spectrum so as to align with the rectifying devices, e.g., rectenna.^{15} The rectenna, a bowtie antenna with an MIM tunneling diode at its feedpoint, utilizes the ordered electromagnetic radiation from the metamaterial and converts this AC induced from the electromagnetic plasmonic waves to an electrical DC.

Our metamaterial consists of copper and a photoresist, SU-8, in alternating micrometer-scale patterns that are fabricated using various techniques, including photolithography, wet and dry etching of metals and oxides, thermal evaporation of thin-layer metals, and electrodeposition of bulk copper. Our previous paper^{16} has demonstrated the design and fabrication of the 3D metamaterial. This study focuses on how metamaterials respond to EM waves in their designed operation range and achieves the optimization of metamaterial geometry. As is known, there are three types of interactions between incident radiation and irradiated object: transmission ($T$), reflection ($R$), and absorption ($A$). The summation of those three portions is unity,

Since absorption is difficult to directly quantify, it can be calculated after quantifying the transmission and reflection. All transmission and reflection simulations are carried out with FEM based on Maxwell’s equation.

## II. THEORY/CALCULATION

A finite element method is a numerical method for solving problems by subdividing the domain of simulation structure into disjoint components of simple geometry called finite elements.^{17} Some frequency-dependent material properties are needed for FEM, including relative permittivity and refractive index. The refractive index is typically expressed as a complex quantity,

where $n~$ is the complex refractive index, $n$ is the real part of refractive index, and $k$, the extinction coefficient, represents the imaginary part of refractive index. Also, the complex relative permittivity is

where $\epsilon r~$ is the complex relative permittivity and $\epsilon R$ and $\epsilon i$ are the real part and the imaginary part of the complex relative permittivity.

Since the relationship between the complex relative permittivity $\epsilon r~$ and the complex refractive index $n~$ is represented as

thus,

In order to obtain $\epsilon R$ and $\epsilon i$, Ordal *et al.*^{18} introduced a Drude model parametrization for transition metals in the near IR. The Drude model dielectric function is represented as

where $\epsilon \u221e$ is the high frequency permittivity contribution that is usually superior or equal to one. $\omega p$ is the plasma frequency and $\omega \tau $ is the damping frequency. Here, $\omega $, $\omega p$, $\omega \tau $ have units of $cm\u22121$. The real part and imaginary part of the complex dielectric function are

The result of a Drude model fit of the metals is shown in Table I. After substituting $\omega \tau $ and $\omega p$, the expression of $\epsilon R$ versus frequency and $\epsilon i$ versus frequency can be obtained. The refractive index $n$ and the extinction coefficient $k$ can then be calculated through the derivation shown in Eqs. (5)–(7).

Metal . | 10^{−2}ω_{τ} (cm^{−1})
. | 10^{−4}ω_{ρ} (cm^{−1})
. |
---|---|---|

Cu | 0.732 | 5.96 |

Al | 6.60 | 11.9 |

Au | 2.15 | 7.28 |

Ni | 3.52 | 3.94 |

Metal . | 10^{−2}ω_{τ} (cm^{−1})
. | 10^{−4}ω_{ρ} (cm^{−1})
. |
---|---|---|

Cu | 0.732 | 5.96 |

Al | 6.60 | 11.9 |

Au | 2.15 | 7.28 |

Ni | 3.52 | 3.94 |

In this study, plane waves are selected as a radiation source to simplify simulations. As is known, the electric field strength of wave is dependent on space and time, which can be written in the form

Thus,

where bold-faced $E$ is the complete electric field, the lighter type $E$ is the time-independent electric field, and $\omega $ is the angular frequency of the wave. The function $E(x,y,z)$ accounts for the spatial dependence of the field, while $ei\omega t$ represents the time variation. Since the vector wave equation^{19} is expressed as

where $\mu 0$ is the vacuum permeability, $\epsilon 0$ is the vacuum permittivity, and $\epsilon r$ is the relative permittivity. So, the time-independent vector wave equation is

For plane waves (only Z direction),

Guess the solution of Eq. (15) is

Because

where $c\xb0$ is the speed of light in vacuum. In addition, the relationship between relative permittivity and refractive index is shown in Eqs. (4) and (5), thus

The time-independent wave equation of plane wave propagating in the $z$ direction is shown in Eq. (20),

Since the refractive index $n$ and the extinction coefficient $k$ are known, the wave equation of plane wave as shown in Eq. (20) can be solved. The spectroscopic transmission and reflection can be represented using $S$-parameters^{20} and the square of the absolute value of $S$-parameter as shown in Eqs. (21) and (23),

Since

where $E0$ is the electric field of incident radiation, $E1$ is the electric field of reflected portion of the radiation, $E2$ is the electric field of transmitted portion of the radiation, $S11$ is the $S$-parameter of reflection, $S21$ is the $S$-parameter of transmission, and $P$ is the electric power of the radiation.

## III. EXPERIMENT

### A. Setup of transmission simulation

The geometry of the 3D transmission simulation unit is shown in Fig. 1. From the top to the bottom, the blocks are the perfect matched layer (PML), port 1, air, metamaterial unit, and port 2. The PML is a perfect absorber that will absorb unwanted reflections in the system. Each metamaterial unit consists of three parts: a silicon substrate, a SU-8 pillar, and a layer of copper. Figure 1(a) shows an isometric projection displaying the general finite element mesh. In Fig. 1(b), a $y$–$z$ planar cross-sectional plane near the center of the mesh displays each of the distinct regions of layered functional metamaterial structure required to simulate the optoelectronic response to terahertz irradiation. For the simulation, a 1 W plane wave propagates from port 1 to port 2. In Fig. 1(b), $d$ is the pillar height, $Ap$ is the bottom side length, $P$ is the spacing between pillars, $d2$ is the thickness of overplated copper, and $d3$ is the distance of the air layer.

The initial parameters of the geometry and the range of the radiation are shown in Table II. The radiation frequency range is from 0.4 to 2 THz. Since there is a requirement of MMs, the feature size must be smaller relative to the external wave we are interested in. Based on our target radiation range whose peak is at 1 THz, i.e., 300 $\mu $m wavelength, all of the feature sizes must be smaller than 300 $\mu $m. The customized mesh size was set relative to the minimum and maximum wavelength of the terahertz radiation.

Name . | Description . | Expression . |
---|---|---|

Ap | Bottom side length | 50 × 10^{−6} m |

P | Pillar spacing | 220 × 10^{−6} m |

d | Pillar height | 150 × 10^{−6} m |

d_{2} | Copper thickness | 40 × 10^{−6} m |

d_{3} | Air thickness | c_const/f_{min}/4 |

f_{min} | Minimum frequency | 0.4 × 10^{12} Hz |

f_{max} | Maximum frequency | 2 × 10^{12} Hz |

fp_{Cu} | Drude model parameter 1 | 1.79 × 10^{15} Hz |

ft_{Cu} | Drude model parameter 2 | 2.19 × 10^{12} Hz |

λ_{min} | Minimum wavelength | c_const/f_{max} |

λ_{max} | Maximum wavelength | c_const/f_{min} |

df | Step size of frequency | 0.1 × 10^{12} Hz |

Si_{thick} | Silicon thickness | 40 × 10^{−6} m |

Name . | Description . | Expression . |
---|---|---|

Ap | Bottom side length | 50 × 10^{−6} m |

P | Pillar spacing | 220 × 10^{−6} m |

d | Pillar height | 150 × 10^{−6} m |

d_{2} | Copper thickness | 40 × 10^{−6} m |

d_{3} | Air thickness | c_const/f_{min}/4 |

f_{min} | Minimum frequency | 0.4 × 10^{12} Hz |

f_{max} | Maximum frequency | 2 × 10^{12} Hz |

fp_{Cu} | Drude model parameter 1 | 1.79 × 10^{15} Hz |

ft_{Cu} | Drude model parameter 2 | 2.19 × 10^{12} Hz |

λ_{min} | Minimum wavelength | c_const/f_{max} |

λ_{max} | Maximum wavelength | c_const/f_{min} |

df | Step size of frequency | 0.1 × 10^{12} Hz |

Si_{thick} | Silicon thickness | 40 × 10^{−6} m |

Relative permittivity ($\mu r$) and relative permeability($\epsilon r$) are the required material properties for the simulation. For air and silicon, the material properties can be found in the COMSOL library. Since the relative permittivity of SU-8 and copper is frequency dependent, their relative permittivities can be extrapolated from the data based on the Drude model,^{18,21} as shown in Table III.

### B. Setup of transmission measurement

In addition to transmission simulation, a series of transmission experiments for metamaterials were carried out with a terahertz time-domain spectrometer in transmission mode. The setup of the THz-TDS system is demonstrated in the previous paper.^{12} The radiation generated from the system ranges from 0.4 to 5 THz. For each transmission measurement, the electric field magnitude versus time domain spectrum is recorded, then the frequency spectrum is obtained by Fourier transform.

### C. Setup of reflection simulation

In Fig. 2, the geometry setup of reflection is the same as that of the transmission simulation except the metamaterial unit is upside down so as to better display the localized surface plasmon at the metal/dielectric interface. There is still a 1 W plane wave propagating from port 1 to port 2.

## IV. RESULTS AND DISCUSSION

### A. Results and discussion of transmission measurement

Figures 3 and 4 are the transmission spectra to metamaterials with different thicknesses of copper grown on the substrates with 150 $\mu $m-tall pillars and 40 $\mu $m-tall pillars, respectively. Comparing Figs. 4 and 3 shows a series of sharper resonance peaks, which indicates that compared to 40 $\mu $m pillars, 150 $\mu $m pillars provide enhanced resonance to MMs at 1 THz.

Also, as shown in Fig. 3, the transmission intensity decreases with an increase in the copper thickness, and a surface plasmon resonance appears around 1 THz. For the metamaterial with 100 nm thick copper, no metamaterial shift was observed at 0.5 THz, where surface plasmon resonance takes place. As the copper thickness increases from 100 nm to 150 $\mu $m, the surface plasmon appears at a higher radiation frequency, from 0.5 to 1.2 THz, while the maximum transmission drops from 100% to 17%. The appearance of transmission peaks indicates that the geometry of metamaterials, especially the pillar structure, has a significant correlation to the surface plasmon resonance.

Besides, in Fig. 3, the transmission completely disappears when the copper layer is 200 $\mu $m thick. As is common knowledge, irradiated objects can absorb, reflect, or transmit the radiation in some combination. If the transmitted portion is zero, the quantity absorbed can be easily obtained by just measuring the reflection spectrum, which simplifies the optical characterization and the geometry design of metamaterials. Comparing all of the different geometric parameters, there are three potential reasons for zero transmission: overplating, metamaterial thickness, and the metal/dielectric height ratio. In order to investigate the cause of zero transmission, two overplated metamaterials were manufactured to compare their tramsmission spectra as shown in Fig. 5. It is clear that the overplated metamaterials with 40 $\mu $m thick pillar do not give rise to zero transmission, indicating that overplating is not the cause of zero transmission. Figures 6–8 show the transmission spectra of MMs with different metal/dielectric height ratios ranging from 0 to 1. It is apparent that all the metamaterials can transmit terahertz radiation, indicating that the metal/dielectric height ratio is not the reason for zero transmission. After ruling out the cause of overplating and metal/dielectric height ratio, we can conclude from Figs. 5–8 that the transmission has a strong correlation with the thickness of metamaterials. In our study, this transmission disappears when the metamaterial consists of 150 $\mu $m thick pillars and 200 $\mu $m thick copper.

As shown in Fig. 6, the transmission performance of MMs with 100 nm thick copper and different pillar heights is very similar in that both of them have a significant transmission spectrum at around 1 THz and that there appears 100% transmission at 0.7 THz. It means that metamaterial pillars have the ability to manipulate terahertz radiation.

As shown in Fig. 7, the metamaterial with both 150 $\mu $m pillars and 75 $\mu $m copper gives rise to a sharp transmission spectrum, showing a surface plasmon resonance at 1 THz with a bandwidth of $\xb10.1$ THz. However, a double transmission peak at $1\xb10.4$ THz appears when the terahertz wave irradiates on the metamaterial with 40 $\mu $m pillars and 20 $\mu $m copper. This indicates that, compared to the metamaterial with 40 $\mu $m pillars and 20 $\mu $m copper, the metamaterial with 150 $\mu $m pillars and 75 $\mu $m copper can couple terahertz radiation to the propagating plasmon and generate surface plasmon resonance at 1 THz more efficiently.

Figure 8 shows the comparison of transmission for metamaterials with the same metal/dielectric heights. In general, the metamaterial with 150 $\mu $m pillar and copper has a weaker transmission spectrum and less noise when compared to metamaterial with 40 $\mu $m pillars and copper. The maximum transmission of metamaterial with 150 $\mu $m thickness is 16%, whereas this is 52% for metamaterial with 40 $\mu $m thickness. Figure 8 indicates that with the same metal/dielectric height, thicker metamaterials give rise to a more coherent surface plasmon and less noise. Also, it proves again that thicker metamaterials cause a decrease in the transmission magnitude.

### B. Results and discussion of transmission simulation

Figures 9 and 10 are the experimental and simulation results of a terahertz transmission test on metamaterials with 150 $\mu $m thick pillars and different thicknesses of copper. In general, the two results show the same trend that the transmission decreases with an increase in the copper thickness, and there appears to be a surface plasmon resonance around 1.2 THz. For the metamaterial with 150 $\mu $m thick pillars and 40 $\mu $m thick copper shown in Fig. 9, a surface plasmon localizes at 1.17 THz with 96% terahertz radiation being transmitted through the metamaterial. As the copper thickness increases to 75 $\mu $m and further increases to 150 $\mu $m, the transmitted radiation peak drops to 90% and 45%, respectively. When the copper thickness increases to 200 $\mu $m, terahertz radiation no longer transmits through the metamaterial. The disappearance of transmission is also observed in the simulation result in Fig. 10. However, there is a significant difference between the experimental and simulation results in that the transmission amplitude in the simulation is lower than that measured in the experiment. This is due to defect-free lattices and perfect lamination between metal/dielectric layers in the FEM model, which gives rise to more reflection and absorption, resulting in less transmission. Moreover, the reflected portion is absorbed by the PML in the FEM model, preventing it from reflecting back. As a result, the amount of the radiation striking the metal/dielectric in the simulation is less than that of the experimental result, which gives rise to a generally weaker transmission spectrum from the simulation.

### C. Results and discussion of reflection simulation

In order to optimize the geometry of the metamaterial, all the parameters in the system are investigated including the pillar height ($d$), bottom side length ($Ap$), spacing between pillars ($P$), the thickness of overplated copper ($d2$), the distance of the air layer ($d3$), and the metal type. As shown in Figs. 11–13, even as we vary the metal type (copper, aluminum, gold, and nickel), the thickness of the air layer (from 190 $\mu $m to 750 $\mu $m), and the thickness of overplated copper (from 0.1 $\mu $m to 40 $\mu $m), there is no change in the terahertz transmission spectrum. All of the spectra show three minimum points: 70% of the radiation is reflected at 1.2 and 1.4 THz, reflectively, and 67% of the radiation is reflected at 1.8 THz. The independence of reflection to these parameters means that metal only acts as a medium to transfer the propagating surface plasmon, so the difference in the electrical conductivity between different metals will not affect the reflection property of the metamaterial.

However, the reflection property of the metamaterial changes significantly if the SU-8-related parameters ($d$, $Ap$, and $P$) are changed, as shown in Figs. 14–16. As shown in Fig. 14, with an increase in the pillar height ($d$) from 50 to 125 $\mu $m, the minimum reflection of the metamaterial shifts from the high frequency to low frequency, indicating that the surface plasmon resonance is also moving from high frequency to low frequency; and the minimum reflection is 39% localized at 1.2 THz when the pillar height is 100 $\mu $m. As shown in Fig. 15, the increase in the bottom side length ($Ap$) from 40 to 70 $\mu $m gives rise to a decrease in terahertz pulse reflection from MMs, and there is a minimum reflection, 63%, taking place at 1.3 THz. In addition, the increase in the spacing ($P$) between pillars from 120 to 200 $\mu $m increases the reflection on metamaterials, as shown in Fig. 16. When we combine the results from Figs. 15 and 16, we found that with the increase in the surface area ratio of SU-8/Cu, the terahertz reflection to MMs decreases.

After considering each of the parameters from all perspectives, including the absorption efficiency of metamaterials, the ease of fabrication and characterization, and the compatibility between metamaterials and the diode, the geometry of our three-dimensional metamaterial is optimized as shown in Fig. 17. The pillar height is 150 $\mu $m, the aperture dimension is 80 $\mu $m, the spacing between pillars is 220 $\mu $m, and the thickness of overplated copper is 50 $\mu $m. In this geometry, the transmission portion is zero, so it is straightforward to characterize the reflection and absorption property of MMs. Copper is the metallic material of choice in the metamaterial due to its high thermal and electrical conductivity, excellent mechanical properties, relative cost, and availability.

Figure 18 shows the reflection spectrum of the MMs with the optimized geometry. It is clear that there are three major surface plasmon resonance peaks observed at frequencies of 1.0 THz (70%), 1.2 THz (70%), and 1.8 THz (68%). The electric field distribution of the optimized metamaterial at different frequencies is shown in Fig. 19. Color mapping is used to display electric field intensity, with red indicating the highest amplitude. The red color represents the highest electric field region. Figure 19(b) shows that the highest electric field is located inside the dielectric pillar and also at the metal/dielectric/air interface, which is also the location to position the nano-antenna in the future; however, the radiation of higher or lower frequency cannot be focused on the MMs. So we conclude that the optimized geometry can localize the surface plasmon with the peak frequency of 1.2 THz and reflection intensity of 70%. This means that the highest terahertz absorption is 30%, which takes place at 1.2 THz.

## V. CONCLUSIONS

A series of FEM COMSOL simulations were carried out to investigate the transmission and reflection properties of metamaterials. The transmission experiments show that the geometry of metamaterial, especially the pillars structure, has a significant correlation to surface plasmon resonance. Compared to the metamaterial with 40 $\mu $m pillars and 20 $\mu $m copper, 150 $\mu $m pillars and 75 $\mu $m copper can couple THz radiation to the propagating plasmon and generate surface plasmon resonance at 1 THz more efficiently. In addition, the simulation and experimental results show the same trend that the transmission decreases with an increase in the copper thickness. Moreover, no terahertz radiation transmits through the metamaterial with 150 $\mu $m thick pillar and 200 $\mu $m thick copper.

The reflection simulation shows that the reflection to metamaterial is dependent on SU-8 related parameters, including pillar height, bottom side length, and spacing between pillars, while the reflection is independent of the thickness of overplated copper, the distance of the air layer, and the metal type. We also found that with the increase in the surface area ratio of SU-8/Cu, the terahertz reflection to MMs decreases. Finally, the optimized metamaterial geometry is achieved and it can localize the surface plasmon with the peak frequency of 1.2 THz to the nano-antenna and the highest terahertz absorption is 30%.

## ACKNOWLEDGMENTS

This work was funded by RedWave Energy, Inc. and supported through a grant from the U.S. DOE-ARPA-E (Grant No. 1261-3407). The authors also wish to thank the support from the MU engineering library.

## DATA AVAILABILITY

The raw data of this research were generated by COMSOL and THz-TDS system at the University of Missouri. The data that support the findings of this study are available from the corresponding author upon reasonable request.