An Eigenmode Study of Nanoantennas from Terahertz to Optical Frequencies

Abstract

In this work, we present a rigorous full-wave eigenanalysis for the study of nanoantennas operating at both terahertz (THz) (0.1–10 THz), and infrared/optical (10–750 THz) frequency spectrums. The key idea behind this effort is to reveal the physical characteristics of nanoantennas such that we can transfer and apply the state-of-the-art antenna design methodologies from microwaves to terahertz and optics. Extensive attention is given to penetration depth in metals to reveal whether the surface currents are sufficient for the correct characterization of nanoantennas, or the involvement of volume currents is needed. As we show with our analysis, the penetration depth constantly reduces until the region of 200 THz; beyond this point, it shoots up, requiring volume currents for the exact characterization of the corresponding radiating structures. The cases of a terahertz rectangular patch antenna and a plasmonic nanoantenna are modeled, showing in each case the need of surface and volume currents, respectively, for the antenna’s efficient characterization.

1. Introduction

Global mobile data traffic grew $71\%$ in 2017, reaching $11.5$ exabytes per month at the end of 2017 [1]. This tremendous growth in mobile traffic will be 77 exabytes per month by 2022; annually, it is believed to reach almost one zettabyte. By 2025, the number of Internet of Things (IoT) worldwide connected devices is expected to reach 75 billion [2]. This immense expansion of wireless devices and wireless data traffic requires extremely low latencies in the order of 1 ms, real-time data transfer, strong determinism and low-delay jitter in the order of several $\mathsf{\mu }$s [3]. Current communication systems operating at and below 60 GHz are not able to satisfy these high needs [4]. Therefore, the investigation of devices operating at higher frequency bands that can satisfy the need for high data rates, allowing massive connectivity between heterogeneous devices, is needed. To this end, the terahertz (THz) frequency band in the range of 0.1–10 THz, and the infrared and optical spectrum, in the range of 10–750 THz, have recently received noticeable attention in the research community [5,6,7,8,9], appearing as ideal choices for high-speed communication links.

Namely, in 2008, the IEEE $802.15$ suite of standards defined regulatory frameworks over the frequency band ranging between 275 GHz and 3000 GHz [10], which was revised to 802.15.3 in 2014, with the goal of corresponding communication systems being able to achieve data rates in the order of 100 GBps. As also listed by the International Telecommunication Union, the main characteristics of terahertz are as follows [11]: (i) safety—non-ionizing photon energy, (ii) high spectral resolution, (iii) high spatial resolution, (iv) high penetration power, and (v) rapid attenuation in water (the last two characteristics are of less importance, unless more emphasis is put on security in wireless communications). Therefore, these characteristics, together with the much larger bandwidth provided at terahertz frequencies, can satisfy the requirements posed by the exponential growth of wireless data traffic. Terahertz links may offer data rates in the order of 1 TBps, in contrast to mmWave systems that offer links in the order of 10 GBps [4]. Furthermore, terahertz signal propagation offers higher directionality and enables secure communications over highly sensitive links, due to the demand of extremely high gain antennas.

Extensive research has been performed over the last decade on mmWave (e.g., [12,13]), terahertz (e.g., [14,15,16,17,18,19,20,21,22,23]), plasmonic (e.g., [24,25,26,27,28,29]) and optical (e.g., [30,31,32,33,34]) antennas, with the researchers being focused on the following: (a) applying established microwave frequencies techniques at these high frequencies (e.g., [35,36]); (b) identifying suitable materials that are appropriate for design at frequencies over 100 GHz (e.g., [18,23]); or (c) identifying appropriate and low-cost fabrication approaches (e.g., [32,37,38,39]). In this work, we are focused on the first thrust.

Even though microwaves techniques are scalable, the transition from RF-microwaves to optical frequencies needs to be done thoroughly. One of the main aspects that needs to be taken into account is the electromagnetic response of metals at high frequencies. At photonic frequencies, metals behave like lossy dielectrics, and they can be described with the introduction of a complex permittivity [40], whose accurate characterization is challenging due to its frequency dependence. Additionally, in the ultraviolet regime, metals’ permittivity exhibits a relaxation from negative to positive values, while at their plasma resonant frequency of about $2.4\times {10}^{15}$ Hz, it becomes transparent to electromagnetic waves [41]. It is due to this behavior of metals that the penetration of fields and currents requires a particular examination. This challenge was addressed either with the use of parameterized models [42], or by conducting measurements [43]. For instance, in [44], Novotny generalized the linear scaling rule to optical frequencies. In [45], Greffet et al. determined the impedance of a nanoantenna using electric dipole moments. Alù and Engheta et al., with a series of works, described and modeled the features of optical nanoantennas [32,46], tuned their scattering responses through loading [47], and utilized the concept of optical circuit elements [48]. Finally, in [49], Zheng et al. established a natural mode analysis for nanotopologies operating at photonic frequencies.

In this work, we focus on the eigenanalysis of nanoantennas for the frequency range 0.1–750 THz, exploiting the offered insightful view of the studied structures. Utilizing the attained eigensolutions (eigenvalues and corresponding eigenvectors), we are able to determine an in-depth characterization of the nanoantennas’ behavior. Namely, the eigenvalues inform us of the frequency of operation of these antennas, while the eigenvectors identify the modal field and current distributions which indicate the appropriate excitation type. Most importantly, the eigenvectors reveal the behavior of the corresponding materials at these frequencies, and in particular, they inform us about the penetration depth in metals, which serves as the objective of this work. The latter is of critical importance for both theoretical and practical considerations. Namely, from the theoretical point of view, it reveals whether surface currents are sufficient for the characterization of nanoantennas, or the involvement of volume currents is needed. From the practical point of view, it suggests the appropriate thickness of metal needed to maximize the efficiency of structures designed at these frequencies. In this study, we particularly show how the efficiency of a nanoantenna varies, as we take into consideration the appropriate metal thickness. Specifically by utilizing the current flowing in the metallic radiating areas, either surface or volume, we observe that currents penetrating in the metallic regions have a small effect on the efficiency of nanoantennas, suggesting in practical applications the use of thin metallic layers.

The proposed formulation is based on our prior developed eigenvalue technique [50], which is appropriate for the study of any arbitrarily shaped geometry loaded with anisotropic and/or inhomogeneous materials. In what follows, the proposed formulation is demonstrated in Section 2; the penetration depth in metals is investigated in Section 3; and Section 4 cross-examines the results of published studies and those derived from a commercial simulator. Section 5 concludes our work, giving future research directions.

2. Formulation

In this work, the finite element method (FEM) is utilized for the eigenanalysis of the studied geometries. Let us assume the source-free, and in general, unbounded $\mathsf{\Omega }$ domain shown in Figure 1. Absorbing boundary conditions (ABCs) of the first kind are used to apply the radiation condition instead of perfectly matched layers (PMLs), avoiding the well-known spurious Berenger modes [51]. Following our previous work in [52], the resulting formulation is general and able to describe any arbitrarily shaped $3D$ computational domain, inhomogeneously loaded with, in general, anisotropic material of tensor permittivity (${\overline{\overline{ϵ}}}_r$), tensor permeability (${\overline{\overline{\mu }}}_r$), and tensor conductivity ($\overline{\overline{\sigma }}$). Perfect electric conductors (PEC), $S_E$, perfect magnetic conductors (PMC), $S_H$, and the open-radiating fictitious surface, $S_{fi}$, where the ABCs are applied, comprise the domain boundary $\mathsf{\partial }\mathsf{\Omega }$ (namely, $\mathsf{\partial }\mathsf{\Omega }=S_E\bigcup S_H\bigcup S_{fi}$). Although the utilized formulation is already given in our previous work [50], let us summarize this approach in the next section for a convenient and self-sustained paper.

2.1. Eigenvalue Problem

Following [53], Maxwell’s differential equations along with the involved boundary conditions, for both the electric and magnetic field, read as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{cc}\overline{\nabla }\times \overline{E}=-j\omega \mu {}_0\overline{\overline{\mu }}{}_r\overline{H}& \overline{\nabla }·(\overline{\overline{ϵ}}{}_r-j\frac{\overline{\overline{\sigma }}}{\omega {ϵ}_0})\overline{E}=0\\ \overline{\nabla }\times \overline{H}=j\omega ϵ{}_0\overline{\overline{ϵ}}{}_r\overline{E}& \overline{\nabla }·\left(\overline{\overline{\mu }}{}_r\overline{H}\right)=0\end{array}\mathrm{in}\mathsf{\Omega }\\ \hfill \widehat{n}\times \overline{E}\times \widehat{n}{|}_{S_E}=0\mathrm{on}\mathrm{PEC}\\ \hfill \widehat{n}\times \overline{H}\times \widehat{n}{|}_{S_H}=0\mathrm{on}\mathrm{PMC}\\ \hfill \widehat{n}\times (\overline{\nabla }\times \overline{E}){{|}_{S_{fi}}\approx -j{k}_0\widehat{n}\times (\widehat{n}\times \overline{E})|}_{S_{fi}}\mathrm{on}\mathrm{the}\mathrm{fictitious}\mathrm{surface}\end{array}\end{array}$$

(1)

where ${ϵ}_0$ and ${\mu }_0$ are the free space permittivity and permeability, respectively. $\omega$ is the operating resonant frequency, which constitutes the eigenvalue to be sought (i.e., resonant frequency); $\widehat{n}$ is the outward unit normal vector defined on the fictitious surface $S_{f_i}$; and $k_0$ is the resonant wavenumber of free space, which is connected with the circular resonant frequency through the relation $k_0=\omega /c$, where c is the speed of light.

For sake of convenience, a complex permittivity accounting for both displacement current and conductivity current losses is introduced as follows:

$${\overline{\overline{ϵ}}}^{\ast }=\overline{\overline{ϵ}}-j\frac{\overline{\overline{\sigma }}}{\omega }={ϵ}_0{\overline{\overline{ϵ}}}_r^{\ast },{\overline{\overline{ϵ}}}_r^{\ast }={\overline{\overline{ϵ}}}_r-j\frac{\overline{\overline{\sigma }}}{\omega {ϵ}_0}$$

(2)

For isotropic media, two alternative representations may be considered. In the first one, all losses are included in a complex dielectric constant ${ϵ}_r^{\ast }={ϵ}_{rR}-j{ϵ}_{rI}$ as in (2), while in the second one, the total losses are accounted within an equivalent conductivity ${\sigma }_{eff}$ as follows:

$$\overline{\nabla }\times \overline{H}=j\omega {ϵ}_0{ϵ}_{rR}\overline{E}+\omega {ϵ}_0{ϵ}_{rI}\overline{E}=j\omega {ϵ}_0{ϵ}_{rR}\overline{E}+{\sigma }_{eff}\overline{E}$$

(3)

Thus, it is ${\sigma }_{eff}=\omega {ϵ}_0{ϵ}_{rI}$.

It should be noted that these types of losses retain the same form of eigenproblem as the lossless case, apart from the fact that the resulting eigenvalues $k_0^2$ become complex. On the contrary, the incorporation of ABCs introduces an additional term $k_0$, besides $k_0^2$, rendering the eigenproblem of a polynomial nature.

Eliminating the magnetic field $\overline{H}$ from Maxwell’s equations in (1), the electric field vector wave equation is formed as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \overline{\nabla }\times \left[\overline{\overline{\mu }}{{}_r}^{-1}·\overline{\nabla }\times \overline{E}\right] - k_0^2\overline{\overline{ϵ}}{}_r\overline{E}=0\end{array}\end{array}$$

(4)

Following the standard Galerkin procedure, (4) is written in its weak form as follows [53]:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\iiint }_{\mathsf{\Omega }}(\overline{\nabla }\times \overline{T})·\overline{\overline{\mu }}{{}_r}^{-1}(\overline{\nabla }\times \overline{E})d\mathsf{\Omega }-k_0^2{\iiint }_{\mathsf{\Omega }}\overline{T}·\left(\overline{\overline{ϵ}}{}_r\overline{E}\right)d\mathsf{\Omega }-j{k}_0{Z}_0{∯}_{Si}\overline{T}·(\widehat{n}\times \overline{H})d{S}_i=0\end{array}\end{array}$$

(5)

where $\overline{T}$ is the vector weighting function, $Z_0=120\pi$ is the free space impedance, and i in $S_i$ can be any or all of $E,H,f_i$ as appears in (1). The surface integral is defined over the surface $dS$ of any element. Applying the 1st kind of ABCs [54], (5) takes the following form:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\iiint }_{\mathsf{\Omega }}(\overline{\nabla }\times \overline{T})·\overline{\overline{\mu }}{{}_r}^{-1}(\overline{\nabla }\times \overline{E})d\mathsf{\Omega }-k_0^2{\iiint }_{\mathsf{\Omega }}\overline{T}·\left(\overline{\overline{ϵ}}{}_r\overline{E}\right)d\mathsf{\Omega }+j{k}_0{∯}_{S_{fi}}(\widehat{n}\times \overline{T})·(\widehat{n}\times \overline{E})d{S}_{fi}=0\end{array}\end{array}$$

(6)

In matrix form, (6) can be written as follows:

$$\begin{array}{c}\hfill \left[S\right]\left[E\right]-k_0^2\left[M\right]\left[E\right]+k_0\left[R\right]\left[E\right]=0\end{array}$$

(7)

with

$$\begin{array}{c}\hfill \left[S\right]={\iiint }_{\mathsf{\Omega }}(\overline{\nabla }\times \overline{T})·{\overline{\overline{\mu }}}_r^{-1}(\overline{\nabla }\times \overline{E})d\mathsf{\Omega }\end{array}$$

(8)

$$\begin{array}{c}\hfill \left[M\right]={\iiint }_{\mathsf{\Omega }}\overline{T}·\left({\overline{\overline{ϵ}}}_r\overline{E}\right)d\mathsf{\Omega }\end{array}$$

(9)

$$\begin{array}{c}\hfill \left[R\right]=j{∯}_{S_{fi}}(\widehat{n}\times \overline{T})·(\widehat{n}\times \overline{E})d{S}_{fi}\end{array}$$

(10)

Overall, (7) is general and can be used for the eigenanalysis of any open radiating structure, even for designs operating at photonic frequencies, where the infinite conductivity approach of metals is no more valid. In this case, metals are modeled as frequency-dependent lossy dielectrics with a negative real part of the dielectric constant, making the matrices of (7) also frequency dependent. Therefore, the polynomial nature of the traditional low microwave frequencies system of (7) transforms into a purely non-linear system, increasing its solution complexity. To address this problem, a linearization approach is utilized as discussed in Section 2.2.

2.2. Photonic System Eigenproblem and Its Linearization Approaches

As mentioned above, when metals are modeled in photonic frequencies, they are mathematically described as frequency-dependent dielectric materials, transforming the polynomial eigenvalue problem of (7) (we assume here that the dielectric materials we use are not frequency dependent, and therefore at low microwave frequencies, (7) is only of polynomial type since S and M have no frequency dependence; note that the wavenumber ($k_0=2\pi f/c$) is the eigenvalue of the problem under solution) to a non-linear eigenvalue problem. To address this non-linearity, a linearization approach is needed. Two different approaches are used: (a) the recursive formulation, and (b) the auxiliary field-based formulation.

2.2.1. Recursive Formulation

In the first approach, an iterative process is introduced, the details of which can be found in our previous work [50]. For reasons of completeness, we only give here a brief description of this method.

We start by setting an initial eigenvalue $k_0^0$ that corresponds to the expected eigenmode. Based on this value, the electric permittivity for all the metals is determined, using the appropriate model (in this work, we use the Lorentz–Drude model [42]), and the mass matrix $\left[M\right]$ is constructed, obtaining an eigenproblem of the form of (20). We solve this problem and we obtain a vector group of complex resonant eigenfrequencies represented as $k_0^1$. An updated mass matrix $\left[M\right]$ is constructed based on $k_0^1$, and a new eigenproblem of (20) form is attained and solved, yielding a new vector group of complex resonant eigenfrequencies, $k_0^2$. The process continues for a number of iterations, $j=1,2,3,\cdots$, obtaining $k_0^j$. It is terminated when two consecutive estimations of the eigenvalue change less than a prescribed error tolerance. Even though this approach is effective, transforming the non-linear eigenproblem to a regular polynomial eigenproblem is very time consuming. Note that the obtained eigenproblem needs to be solved multiple times for each eigenvalue to acquire the desired relative error. Therefore, this linearization approach is efficient only when a small number of eigenvalues is desired.

2.2.2. Auxiliary Field-Based Formulation

When a significantly larger number of eigenvalues is sought, we follow a different approach, recently introduced in both finite difference time domain [55,56] and finite difference [57], finite element [58], and frequency domain simulations for dispersive media. In this second method, a “hidden” auxiliary field is introduced by considering non-magnetic mediums whose permittivities can be modeled by a standard $M–$pole Lorentz–Drude relationship [57,59]:

$$ϵ\left(\omega \right) ={ϵ}_{\infty }\left(1-\sum _{i=1}^M\frac{f_i{\omega }_p^2}{{\omega }^2-{\omega }_{0,i}^2+j\omega {\mathsf{\Gamma }}_i}\right)$$

(11)

In (11), M is the number of oscillators (poles) with frequency ${\omega }_{0,i}$, $f_i$ is the oscillation strength, and ${\mathsf{\Gamma }}_i$ is the relaxation rate, which describes the effective electron scattering rate $({\tau }_i=1/{\mathsf{\Gamma }}_i)$. The parameter ${\omega }_p=(N{e}^2/{ϵ}_{0}m)$ stands for the volume plasma frequency, where N is the electron’s density, e is the electron’s charge, and m is the effective mass of the electron. Finally, the parameter ${ϵ}_{\infty }$ is called instantaneous or high-frequency relative permittivity, and it accounts for the net contribution from the positive ion cores. The typical values for ${ϵ}_{\infty }$ for the utilized metals range between 1 and 10, depending on the intraband response. For ideal free electron gas, it takes the value ${ϵ}_{\infty }=1$, which is a common choice in many published works.

Let us consider the case of a single-pole Lorentz–Drude model, and introduce two auxiliary fields as follows [59]:

$$\overline{P}=-{ϵ}_{\infty }\frac{f_0{\omega }_p^2}{{\omega }^2-{\omega }_0^2+j\omega \mathsf{\Gamma }}\overline{E}$$

(12a)

$$\overline{J}=j\omega \overline{P}$$

(12b)

where $\overline{P}$ is the polarization field, which is correlated with the electron’s concentration as $\overline{P}=Ne\overline{r}$ ($\overline{r}$: position), and $\overline{J}$ is the current density. Thus, the basic equations of motion for the electromagnetic field in a dispersive (non-magnetic) medium are as follows [57]:

$$\frac{\mathsf{\partial }\overline{H}}{\mathsf{\partial }t}=-\frac{1}{{\mu }_0}\overline{\nabla }\times \overline{E}$$

(13a)

$$\frac{\mathsf{\partial }\overline{E}}{\mathsf{\partial }t}=\frac{1}{ϵ}(\overrightarrow{\nabla }\times \overline{H}-\overline{J})$$

(13b)

$$\frac{\mathsf{\partial }\overline{P}}{\mathsf{\partial }t}=\overline{J}$$

(13c)

$$\frac{\mathsf{\partial }\overline{J}}{\mathsf{\partial }t}=f_0{\omega }_p^2{ϵ}_{\infty }\overline{E}-{\omega }_0^2\overline{P}-\mathsf{\Gamma }\overline{J}$$

(13d)

Combining (13a) with (13b) and (13c) with (13d), we obtain a system of coupled wave equations involving the electric field $\overline{E}$ and the electric polarization $\overline{P}$ as follows:

$$\overline{\nabla }\times (\overline{\nabla }\times \overline{E})-k_0^2{ϵ}_{\infty }\overline{E}-k_0^2{ϵ}_0^{-1}\overline{P}=0$$

(14a)

$$k_0^2\overline{P}-j{k}_0\frac{\mathsf{\Gamma }}{c}\overline{P}-\frac{{\omega }_0^2}{c^2}\overline{P}+\frac{f_0{\omega }_p^2{ϵ}_{\infty }}{c^2}\overline{E}$$

(14b)

In matrix notation, (14) reads as follows:

$$\left(\left[\begin{array}{cc}\overline{\nabla }\times \overline{\nabla }\times & 0\\ {ϵ}_{\infty }\frac{f_0{\omega }_p^2}{c^2}& -\frac{{\omega }_0^2}{c^2}\end{array}\right] + k_0\left[\begin{array}{cc}0& 0\\ 0& -j\frac{\mathsf{\Gamma }}{c}\end{array}\right] + k_0^2\left[\begin{array}{cc}-{ϵ}_{\infty }& -{ϵ}_0^{-1}\\ 0& 1\end{array}\right]\right)\left[\begin{array}{c}\overline{E}\\ \overline{P}\end{array}\right] = \left[\begin{array}{c}0\\ 0\end{array}\right]$$

(15)

The system (15) refers to a polynomial eigenvalue problem. The aforementioned approach can be extended for the case of the multiple-pole Lorentz–Drude model by simply increasing the number of auxiliary fields, namely, taking into account more terms from the expansion in (11). Eliminating the frequency of the pole (${\omega }_0=0$), the Lorentz–Drude model (11) reduces to the simple free-electron Drude model. The latter requires the introduction of a single auxiliary field to convert the non-linear eigenproblem into a linear one.

Applying the standard Galerkin procedure and using the 1st kind of ABCs, the weak formulation of system (15) becomes the following [53]:

$$\left(\left[\begin{array}{cc}\left[S\right]& 0\\ \left[P_2\right]& \left[P_3\right]\end{array}\right] + k_0\left[\begin{array}{cc}-j\left[R\right]& 0\\ 0& \left[P_4\right]\end{array}\right] - k_0^2\left[\begin{array}{cc}\left[M\right]& \left[P_1\right]\\ 0& \left[-P_1\right]\end{array}\right]\right)\left[\begin{array}{c}\overline{E}\\ \overline{P}\end{array}\right] = \left[\begin{array}{c}0\\ 0\end{array}\right]$$

(16)

The stiffness $\left[S\right]$, mass $\left[M\right]$, and radiation $\left[R\right]$ matrices are defined in (8)–(10), respectively, while the auxiliary field–related matrices read as follows:

$$\left[P_1\right] ={\iiint }_{\mathsf{\Omega }}(\overline{T}·\overline{P})dV$$

(17a)

$$\left[P_2\right] =\frac{f_0{\omega }_p^2}{c^2}{\iiint }_{\mathsf{\Omega }}(\overline{T}·\overline{E})dV$$

(17b)

$$\left[P_3\right] =\frac{{\omega }_0^2}{c^2}{\iiint }_{\mathsf{\Omega }}(\overline{T}·\overline{P})dV$$

(17c)

$$\left[P_4\right] =-j\frac{\mathsf{\Gamma }}{c}{\iiint }_{\mathsf{\Omega }}(\overline{T}·\overline{P})dV$$

(17d)

The system (16) can be rewritten in a more compact form as follows:

$$\left(\left[M_0\right] + k_0\left[M_1\right] - k_0^2\left[M_2\right]\right)\left[x\right] = \left[0\right]$$

(18)

As it turns out, after the linearization procedure system, (18) has the same form with (7); therefore, both can be solved following the same approach as described in Section 2.3 below.

2.3. Polynomial Linearization

Formulation (7) is of the second-order polynomial form with respect to the eigenvalue $k_0$. Thus, the incorporation of a linearization procedure is needed in order to solve this eigenproblem. A symmetric linearization approach is utilized herein, taking advantage of the symmetric nature of the involved matrices.

By setting $k_0=\lambda$, (7) reads as a single standard quadrature problem with a characteristic polynomial:

$$\begin{array}{c}\hfill P\left(\lambda \right)=-{\lambda }^{2}M+\lambda R+S=0\end{array}$$

(19)

Defining an auxiliary vector $u=\lambda q$, where q refers to the “seeding” vector, and, therefore, using the eigenvector transformation of the form $k_0\left[q\right]=\left[u\right]$, (7) becomes the following [60]:

$$\begin{array}{c}\hfill \left(\lambda \left[\begin{array}{cc}0& -M\\ -M& R\end{array}\right] + \left[\begin{array}{cc}M& 0\\ 0& S\end{array}\right]\right)\left[\begin{array}{c}u\\ q\end{array}\right] = \left[\begin{array}{c}0\\ 0\end{array}\right]\end{array}$$

(20)

As it occurs, (20) is a linear generalized eigenvalue problem, where its eigenvalue corresponds to the complex wavenumber $k_0=\omega /c$, with $\omega ={\omega }_r+j{\omega }_i$. To solve this problem, the Arnoldi shift and invert algorithm is utilized [61]. After the solution of our system, the resonant frequency of the structure is derived, using the real part of the wavenumber as $f_{res}={\omega }_r/\left(2\pi \right)$, while the imaginary part is used to determine the corresponding quality factor as $Q_{res}={\omega }_r/\left(2{\omega }_i\right)$.

3. Investigation of Field Penetration in Metals

A key aspect on the characterization of antennas’ performance at terahertz and optical regimes is the skin depth of the metallic regions. The skin depth depends on the frequency of operation and the conductivity of the material. Its value can be analytically calculated as follows [62]:

$$\begin{array}{c}\hfill {\delta }_s=\frac{1}{\alpha }\end{array}$$

(21)

where $\alpha =Re\left\{\gamma \right\}$, and $\gamma$ is the complex propagation constant of an incident plane wave given as follows:

$$\begin{array}{c}\hfill \gamma =\alpha +j\beta =j\omega \sqrt{{\mu }_0{ϵ}_0{\overline{\overline{ϵ}}}_r^{\ast }}=j\omega \sqrt{{\mu }_0{ϵ}_0{\overline{\overline{ϵ}}}_r}\sqrt{1-j\frac{\overline{\overline{\sigma }}}{\omega {ϵ}_0{\overline{\overline{ϵ}}}_r}}\end{array}$$

(22)

Considering the definition of (3), the skin depth for isotropic metals can be written in terms of either an effective conductivity, or the imaginary part of the complex dielectric constant as follows:

$${\delta }_s=\sqrt{\frac{2}{\omega {\mu }_0{\sigma }_{eff}}}=\sqrt{\frac{2}{k_0^{2}Im\left({ϵ}_r^{\ast }\right)}}$$

(23)

The values of ${\overline{\overline{ϵ}}}_r^{\ast }$ and ${ϵ}_r^{\ast }$ (see (2)) are calculated either using the Drude or single-pole Lorentz–Drude model. Both models give the same distributions in the terahertz frequency spectrum, as the intraband transition of electron is the dominant mechanism, and the contribution of the interband counterpart is negligible [42]. Note that in this study, all modeled structures are non-magnetic; therefore, ${\mu }_r$ is omitted.

As discussed in the introduction, the goal of this work is to investigate whether the surface currents are sufficient for the correct characterization of nanoantennas, or the involvement of volume currents is needed. To do so, we study the penetration depth in conductive regions consisting of copper (Cu), silver (Ag), and gold (Au) for the two frequency ranges of interest: (a) 0.1–10 THz, and (b) 10–750 THz. Note that copper and silver are extensively used at terahertz frequencies [14,15,19,21,22,23], while gold (Au) films are mostly used in photonic and plasmonic structures [24,26,28,31], even though it appears in some works at a lower terahertz, especially in sub-terahertz applications, where gold films are a common choice [20]. The question to be addressed is whether for practically utilized metallic film thicknesses, the electric field penetrates to some significant depth, where the corresponding current density $(\overline{J}=j\omega {ϵ}_0{\overline{\overline{ϵ}}}_r^{\ast }\overline{E})$ must be treated as a volume distribution. Otherwise, if the current epidermal nature prevails, this distribution can be considered to be a surface one (infinitely thin). In the latter case, the rules and approaches established in the microwave and millimeter wave regimes can be directly applied.

3.1. Penetration Depth at 0.1–10 THz

First, we start investigating the penetration depth in the 0.1–10 THz frequency range. As an example, we use a typical patch antenna. The typical thickness (t) of metallic portions for a patch antenna (i.e., the metallic isle on the top and its ground plane) varies from 0.25 $\mathsf{\mu }$m to 0.5 $\mathsf{\mu }$m [21,22], or in some thicker alternatives from 40 $\mathsf{\mu }$m to 50 $\mathsf{\mu }$m [14,15], while known substrate materials like Rogers RO3010, −3006, and −3003 are also made with $t=9-70$ $\mathsf{\mu }$m [63]. (For metal thicknesses below 9 $\mathsf{\mu }$m, lithographic techniques are used (e.g., electron beam lithography (EBL) and focused ion beam (FIB) milling) [64].) The specific thickness is selected based on the maximum (nominal) current density to be encountered (thus, the structure’s operating power) and the specified mechanical rigidity [65]. The penetration depth in metals is inversely proportional with frequency when the metal’s conductivity and dielectric constant vary according to the Drude model (11). Thus, as frequency increases, the penetration depth decreases as depicted in Figure 2. Therefore, the penetration depth is significantly smaller than the typical metal thickness at the high end of the observable frequency range, and the skin effect met in a low-microwave (i.e., $f=1-100$ GHz) patch antenna is still valid at these frequencies. In this case, metals can still be considered perfect electric conductors, and the antenna behavior can be characterized by surface currents flowing on the surface of the patch.

Although surface currents seem to be sufficient for antennas characterization at the 0.1–10 THz frequency range, when the metallic layers are thick, it might not be the case for thicknesses below 2 $\mathsf{\mu }$m. Figure 3, Figure 4 and Figure 5 present a 3D plot of the ratio of the metal thickness $t=0.25-2$ $\mathsf{\mu }$m to the penetration depth ${\delta }_s$, $r=t/{\delta }_s$ at the terahertz frequency range of 0.1–10 THz for the cases of Cu, Ag, and Au, respectively. It is clearly evident in these plots that patch antennas, whose metallic layers thicknesses are in the order of micrometers, exhibit penetration depths of several orders of magnitude lower than the metal thickness. Only for the frequency range below 2 THz, the penetration depth remains comparable to the thickness of the metallic layers. Therefore, for cases of metal thickness $t=0.25-2$ $\mathsf{\mu }$m and for frequencies below 2 THz, volume currents must be used to appropriately characterize the corresponding antennas. A recent report on the design and simulation of terahertz-integrated antennas [66] stated that the increase in conductor thickness results in better radiation efficiency. Thus, in practical situations, thicker metallic claddings are preferred, as was also shown in [67] with a microstrip antenna array operating at 0.1 THz.

3.2. Penetration Depth at 10–750 THz

In our second study, we investigate the penetration depth in metals for the frequency range of 10–750 THz. Only the penetration depth in Ag and Au is studied, as Cu is not met in optical configurations. Figure 6 depicts the penetration depth ${\delta }_s$ in nm and the effective conductivity ${\sigma }_{eff}$ in S/m for these two cases, utilizing both Drude and Lorentz–Drude models as well as the results derived from the measurements of two silver samples [68]. As expected, the penetration depth responses are different for the two models since the contribution of the intraband transition of electron becomes comparable with the interband counterpart [42]. However, the penetration depth of Ag based on measurements has a response closer to the simple Drude model.

Regarding the thickness of the metallic portions for optical nanoantennas, it varies at $t=20-100$ nm [24,25,69,70]. (A thickness of $t=20$ nm was used for a gold patch in an experimental setup in [25], while thicknesses for both silver and gold patches at $t=40-100$ nm are found in [24,69,70]). See in Figure 6a how the penetration depth is comparable to the thickness of the nanoantennas’ metallic parts, especially when the interband contribution is taken into account through the Lorentz–Drude model. This is also shown in Figure 7 and Figure 8, which present a 3D plot of the ratio of the metal thickness $t=20-100$ nm to the penetration depth ${\delta }_s$, $r=t/{\delta }_s$, at the frequency range of 10–750 THz. The same situation also exists when the simple Drude model is utilized, due to the small thickness of the metallic portions used at high terahertz frequencies. Therefore, the displacement and polarization currents, flowing along and inside the nanoantenna, become dominant, indicating that the antenna’s characterization should be done through the study of a volume current distribution.

4. Numerical Results

The cases of a terahertz rectangular patch antenna and a plasmonic nanoantenna are modeled in this Results section. As proved next, the surface currents can sufficiently characterize the terahertz patch antenna, while volume currents are needed for the plasmonic nanoantenna. An in depth study is performed for the latter, showing the importance of metal thickness, as it affects the radiation efficiency of the structure.

4.1. Rectangular Patch Antenna at 0.1–3 THz

Let us examine the case of the simple terahertz rectangular patch antenna depicted in Figure 9 with $W_s=147.83$ $\mathsf{\mu }$m, $L_s=119.33$ $\mathsf{\mu }$m, $h=21.5$ $\mathsf{\mu }$m, $W_p=83.33$ $\mathsf{\mu }$m, $L_p=54.83$ $\mathsf{\mu }$m, $t=3$ $\mathsf{\mu }$m, designed on the transparent polyimide substrate material met in terahertz applications [22] with a dielectric constant of ${ϵ}_r=3.5$. (Since we investigate the design and performance of a terahertz patch antenna, the substrate material is assumed to be lossless, and thus, $tan\delta =0$.)

The eigenproblem of the terahertz patch antenna without the feeding line (only the rectangular patch denoted by the red line in Figure 9, the substrate and the ground plane) is solved with the recursive formulation (Section 2.2.1), utilizing a fixed value of silver’s permittivity, determined with the Drude model (11) at 1 THz. The first 10 eigensolutions are obtained and tabulated in the second column of

Table 1. Patch antenna (Figure 9) eigenfrequencies with ($f_{rAL}$) and without ($f_{rR}$, $f_{rA}$) feeding line (red contour) versus those in the absence of patch radiator and feeding line ($f_{rAS}$).

No	${}^1{\mathit{f}}_{\mathit{rR}}$ [THz]	${}^2{\mathit{f}}_{\mathit{rA}}$ [THz]	${}^3{\mathit{f}}_{\mathit{rAS}}$ [THz]	${}^4{\mathit{f}}_{\mathit{rAL}}$ [THz]
1	0.770 + j0.330	0.781 + j0.387	0.799 + j0.433	0.778 + j0.332
2	0.802 + j0.051	0.813 + j0.059	0.815 + j0.058	0.793 + j0.036
3	0.947 + j0.386	0.960 + j0.453	0.958 + j0.440	0.951 + j0.371
4	1.038 + j0.103	1.037 + j0.088	-	1.040 + j0.084
5	1.171 + j0.170	1.187 + j0.200	-	1.181 + j0.200
6	-	-	-	1.314 + j0.042
7	1.408 + j0.089	1.439 + j0.092	1.441 + j0.083	1.438 + j0.109
8	1.506 + j0.101	1.495 + j0.777	1.501 + j0.847	1.490 + j0.098
9	1.586 + j0.691	1.590 + j0.628	1.589 + j0.630	1.571 + j0.725
10	1.617 + j0.731	1.717 + j0.321	-	1.685 + j0.636

${}^1{f}_{rR}$: Eigenfrequencies for the patch antenna of Figure 9 without the feeding line, based on the recursive formulation of Section 2.2.1; ${}^2{f}_{rA}$: Eigenfrequencies for the patch antenna of Figure 9 without the feeding line, based on the auxiliary field-based formulation of Section 2.2.2; ${}^3{f}_{rAS}$: Eigenfrequencies for the finite grounded substrate of Figure 9 in the absence of the patch radiator and the feeding line, based on the auxiliary field-based formulation of Section 2.2.2; ${}^4{f}_{rAL}$: Eigenfrequencies for the patch antenna of Figure 9, based on the auxiliary field-based formulation of Section 2.2.2.

. To avoid any misunderstanding, recall that eigenmodes are independent of any excitation. However, their modal current (magnetic field) or voltage (electric field) distribution may identify the appropriate position of the excitation along with which the quantity, current or voltage can be maximized at the corresponding resonance. In turn, these resonances are characterized as parallel (maximum voltage at the feed) and series (maximum current at the feed) ones.

The eigenvalues of the terahertz patch antenna are cross examined with those (third column in Table 1) resulting from the eigenanalysis of the same structure, but utilizing the auxiliary field-based formulation (Section 2.2.2). The small deviations of the average order of $1.3\%$ between the resonant frequencies (real part of eigenvalues), resulting from the two formulations, are due to the different approach used for characterizing the metallic layers. The first approach (recursive formulation in Section 2.2.1) uses a constant permittivity value for the characterization of metals throughout the whole frequency spectrum, while the second one (auxiliary field-based formulation in Section 2.2.2) changes the permittivity of metallic regions according to the attained resonant frequency. It should be noted that CST also uses a fixed permittivity for metals at low terahertz frequencies (below 150 THz) [71].

Closely examining Table 1, one may encounter a plethora of nearby located resonant frequencies. This is unexpected for patch antennas and may lead to the suspicion of pseudo-spurious eigensolutions. However, the experience in printed antenna eigenproblems (especially in their characteristic modes analysis) reveals that the majority of the obtained resonances are due to the finite grounded dielectric substrate. From this point of view, two separate families of eigensolutions are expected: (a) one due to the finite grounded substrate (in the absence of the patch radiator), and (b) one generated by the patch radiator itself (e.g., a patch on infinite substrate). Therefore, the observed eigensolutions for the patch antenna of Figure 9 are the results of mode coupling between these two families. In order to shed light to these phenomena, an eigenanalysis is performed, obtaining the eigensolutions of the structure shown in Figure 9 in the absence of the patch and the feeding line. Therefore, eigenfrequencies $f_{rAS}$ ($1,2,3,7,8,9$) shown in Table 1 are the pure resonances of our finite substrate, which, as expected, are slightly shifted to the corresponding eigenfrequecies when the patch radiator is present (see $f_{rA}$ in Table 1). It is next proved that for these substrate modes, the eigencurrents have strong intensity on the ground plane and are negligible on the patch radiator. On the contrary, eigenfrequencies ($4,5,10$) of the patch radiator are expected to have eigencurrents with strong intensity on the patch itself.

The presence of the excitation, especially the inset feed of Figure 9 modifies the metallic patch geometry and is expected to shift its resonant frequencies. For this purpose, the eigenproblem of our THz patch antenna as depicted in Figure 9 (with the presence of the feeding line) is solved once again, utilizing the auxiliary field–based formulation. The eigenvalues of the latter are listed in Table 1 as $f_{rAL}$. As it can be seen from them, the presence of the feeding line not only causes a subtle change in the eigenfrequencies (compared to $f_{rA}$ in Table 1), but more importantly, it creates a new eigensolution No. $\#6$ as it appears in Table 1. Notably, by carefully observing the patch antenna in Figure 9, the slots of the inset-feed are expected to force the $\widehat{x}–$directed current density ($J_x$) to flow around them, thus increasing their resonant path. Although these modes cannot be excited by this feed, their resonant frequency will be shifted downwards. This is observed, for instance, for modes ($5,10$) of the last column of Table 1.

To prove the validity of the eigenanalysis, Figure 10 presents the real and imaginary parts of both input impedance $Z_{11}$ and input admittance $Y_{11}$ versus the frequency for our terahertz rectangular patch antenna as simulated in CST Microwave Studio [71]. It is by now well established that the true nature of the patch radiator is revealed at its actual input, identified at the feeding line-patch interconnection. This is usually obtained through the de-embedding approach, which is also available in all commercial simulators. Observing the de-embedded input impedance in Figure 10, one may identify multiple resonances at frequencies $1.178,1.705,2.375,2.653$ THz, where $Re\left(Z_{in}\right) \to max$, with an average value of $<Re\left(Z_{in}\right){>}_{av} =139$$\mathsf{\Omega }$. The first two resonances are identified in Table 1 as modes $\#5$ and $\#10$, while the others fall outside the region spanned by the eigenanalysis. Similarly, multiple resonances at frequencies $0.989,1.296,2.033,2.42,2.983$ THz, where $Re\left(Y_{in}\right) \to max$, are observed in Figure 10, with an average value of $<Re\left(Y_{in}\right){>}_{av} =0.0025$ S. Likewise, the first two resonances correspond to modes $\#4$ and $\#6$ of Table 1. All the aforementioned resonances can be identified as different $T{M}_{mn0}$ modes.

It is important to realize that the resonances appearing in the diagram of $Z_{in}$ versus frequency (blue solid line in Figure 10a) are only a subset of the supported modes; namely, these are only modes excited by the specific feed. Additionally, the commonly utilized operating modes in patch antennas are only the parallel resonances appearing at $Re\left(Z_{in}\right)\to max$ (vertical green dotted lines). However, there are always series resonances between successive parallel modes, where $Re\left(Y_{in}\right)\to max$ ($Y_{in}=1/Z_{in}$) (vertical brown dotted lines) as seen in Figure 10a. These series modes are found to be inappropriate for operating–exciting a patch antenna. Hence, the eigenanalysis provides a plethora of antenna modes (see Table 1 besides the substrate modes), which can be discriminated as series and parallel resonances, while only a subset of them can be excited by the feed utilized in Figure 9.

As explained above, it is important to illustrate the modal current distributions in order to identify the position and the type of possible excitations. Surface currents ${\overline{J}}_s$ result from the so-called internal magnetic field (below the patch) as ${\overline{J}}_s=\widehat{n}\times \overline{H}=-\widehat{z}\times \overline{H}$. The eigenvoltage distributions, defined between the patch and the ground plane, can be also illustrated. This can be easily done since $\overline{E}=-\overline{\nabla }V$. Then, the voltage $V(x,y)$ is computed as $V(x,y)=-{\int }_0^h{E}_z(x,y,z)$, which for a thin substrate of height h the ${\overline{E}}_z$ component is nearly homogeneous in the z-direction; thus, $V(x,y)\approx -h{E}_{z}x,y,z$. Nevertheless, since the parallel resonances are only considered operating modes, the series resonances are neglected, which can be simply identified as those whose current distribution at the feed is maximal.

The eigencurrent distribution of mode Nos. 1, 4, 5, and 10—see column $\#3$ in Table 1—flowing upon the top surface of the patch (top isle) and the ground layer are illustrated in Figure 11. The surface current distribution of mode $\#1$ (see Figure 11a) clearly reveals that this is a non-radiating finite substrate mode since its current density is higher on the ground plane rather than on the patch radiator. In turn, the surface current distributions of mode Nos. 4, 5, and 10 correspond to the $T{M}_{010}$, $T{M}_{110}$, and $T{M}_{100}$ radiating patch modes, respectively (see Figure 11b–d). Namely, the current distribution of the dominant $T{M}_{010}$ mode (Figure 11b) is the same with the one observed in the source-driven CST analysis of a typical microwave rectangular patch antenna. Regarding the other (such as No. $\#6$) mode distributions of our patch, they cannot be computed with this source-driven CST analysis since these modes cannot be excited with the specific feeding line. Figure 12 shows the eigencurrent distributions for mode Nos. 3 and 6—column $\#5$ in Table 1—where the modification of distributions due to the presence of the feeding line is clearly evident. Special attention should be given to the current distribution of mode $\#6$ (Figure 12b), which represents a series resonance (maximum current at the feed) and clearly reveals why series resonances are inappropriate for operating–exciting a patch antenna. An insight is offered by Figure 12 regarding the affected current density components. Clearly, the $\widehat{y}–$directed current ($J_y$), which is parallel to the inset feed slots, is negligibly affected, and thus, without any noticeable shift in the resonant frequencies of the corresponding mode Nos. 4, 5, and 10. On the contrary, modes with dominant $\widehat{x}–$directed currents, which are normal to those slots, are strongly affected or even revealed by the corresponding eigenanalysis ($f_{rAL}$ in Table 1).

Overall, it is justified, herein, that the behavior of a patch antenna operating at terahertz frequencies can be once again characterized by the surface currents flowing on its metallic layer.

4.2. Plasmonic Nanoantenna at 400–600 THz

In the second example, we study the plasmonic nanoantenna shown in Figure 13 introduced by Liu et al. [24]. This structure consists of two elliptical cylinder shape gold nanoparticles placed upon a quartz substrate with geometrical characteristics as depicted in Figure 13. The thickness of the particles is $t=40$ nm each and their in-between gap is $g=28$ nm, while their major and minor axes are $r_{max}=108$ nm and $r_{min}=58$ nm, respectively.

Table 2. Ten first eigenvalues and corresponding resonant wavelengths for the plasmonic nanoantenna (see Figure 13) using the auxiliary field-based formulation (Section 2.2.2). The last two columns refer to the eigenvalues and corresponding resonant wavelengths of the substrate (absence of gold particles). The eigenvalues appearing in bold ($3,4$) correspond to patch antenna modes.

	Plasmonic Nanoantenna		Substrate
No	${\mathit{k}}_{\mathbf{0}}$  $\left[{\mathbf{10}}^{\mathbf{6}}\mathit{m}\right]$	$\mathbf{\lambda }$  $\left[\mathit{nm}\right]$	${\mathit{k}}_{\mathbf{0}}$  $\left[{\mathbf{10}}^{\mathbf{6}}\mathit{m}\right]$	$\mathbf{\lambda }$  $\left[\mathit{nm}\right]$
1	9.451 + j4.102	664.8	9.425 + j4.061	666.7
2	9.488 + j3.693	662.2	9.437 + j3.950	665.8
3	9.506 + j3.899	661.3	-	-
4	9.938 + j4.197	632.2	-	-
5	12.210 + j3.273	514.6	12.166 + j3.545	516.5
6	12.226 + j3.403	513.9	12.193 + j3.395	515.3
7	12.290 + j3.649	511.2	12.202 + j3.369	514.9
8	12.321 + 3.656	510.0	12.264 + j3.672	512.3
9	12.395 + j3.310	506.9	12.419 + j3.452	505.9
10	12.455 + j3.428	504.5	12.445 + j3.498	504.9

presents 10 eigenvalues and the corresponding resonant wavelengths of this plasmonic nanoantenna in the frequency range of interest. The same approach as that for the terahertz patch antenna of the previous section is utilized. Hence, the eigenanalysis using the auxiliary field–based formulation (Section 2.2.2) is first carried out for the whole antenna of Figure 13, while this is repeated afterward in the absence of the gold nanoparticles. These eigenvalues are cross examined among them, where it is clearly revealed that only eigenvalues Nos. 3 and 4 are related with the nanoparticles (antenna) modes, while the rest correspond to substrate modes. The electric current density ($\overline{J}$) defined in (12b) includes both the conduction current ($\sigma \overline{E}$) and the displacement current ($j\omega \overline{\overline{ϵ}}\overline{E}$) through the complex permittivity (2) or (11) according to (12b):

$$\overline{J}=j\omega \overline{P}=j\omega {\overline{\overline{ϵ}}}^{\ast }\overline{E}$$

(24)

On the contrary, at lower frequencies, metals behave as perfect electric conductors (PEC), where the tangential electric field vanishes, leading to the assumption of the surface current density ${\overline{J}}_s=\widehat{n}\times \overline{H}$. The resulting current distributions are estimated and plotted in Figure 14.

The two antenna modes ($3,4$) are compared in terms of their eigencurrents and against an indicative substrate mode (1) as illustrated in Figure 14. A close examination of Figure 14c shows that mode No. 4 represents the usual dipole antenna, as its current flows parallel to the long axis of the ellipse and toward the same direction. Conversely, the antenna mode No. 3 exhibits a surprising, strange behavior since its current is normal to the major axis of the ellipse having opposite directions. As these currents are very close, they are expected to yield negligible radiation. Moreover, Figure 14a shows the electric current density flowing upon the gold particles for the substrate mode $\#1$, where its weak intensity, due to the fact that it is a non-radiating mode, is clearly evident. Namely, this substrate mode exhibits a dipole-like current, being parallel to the ellipse’s major axis and flowing in the same direction, but with its intensity being zero at the two edges of the gap. Hence, this mode cannot be excited by the usual, for instance, gap voltage in transmit function, while a detector placed across the gap will receive no response in the receive function. Overall, the eigenanalysis revealed that the useful antenna mode is $\#4$ illustrated in Figure 14c, exhibiting the expected dipolar polarization and resonating at a wavelength of 632 nm, which is in accordance with the original study performed in [24].

Figure 15 shows the real part of the eigenvector at an $x-y$ view 20 nm above the substrate (at the middle-plane crossing the particles) at a wavelength of $632.2$ nm (eigenvalue $\#4$). In this wavelength, the high enhancement of the local field in the gap between the two gold particles, a common characteristic of nanoantenna structures [24], is clearly evident. Additionally, the electric field vector directions reveal the dipolar nature of mode (4).

As stated throughout this paper, a major task of this work is to examine the current distribution over the entire volume of the nanoantenna in order to gain an insight into their radiation mechanism. For this purpose, a side view of the antenna mode distributions for the electric current density is shown in Figure 16. It is clear that the current flows inside the gold nanoparticles, and is not restricted on the surface as happens to antennas operating at low-microwaves (below 100 GHz). Hence, the correct characterization of a nanoantenna operating at optical frequencies needs the use of volume current. Moreover, the current flow is mainly concentrated at the lower interface of the gold layers and the substrate, and is parallel to the substrate for operation at antenna modes (Figure 16b,c). On the contrary, the current flow is concentrated at the top side corner of the gold particles and is normal to the substrate for operation at substrate mode (Figure 16a).

To show the effect of volume currents, in the following subsection, the radiation efficiency of this plasmonic nanoantenna is evaluated as a function of the current’s penetration depth in the gold nanoparticles.

Plasmonic Nanoantenna Radiation Efficiency—Volume Current

As shown in Section 3.2 for the plasmonic nanoantenna of Figure 13, at frequencies of operation above 200 THz, the penetration depth in metals is significant compared to the metal’s thickness. Notably, for the three different cases, at $664.8$ nm, $661.3$ nm and $632.2$ nm, the electric current penetrates in the gold particles, forming the corresponding spatial mode distributions (Figure 16). In order to reveal the effect of the volume currents, the nanoantenna’s radiation efficiency, ${\eta }_R$, is evaluated according to [72] for different thicknesses of metal and for the two antenna modes. Namely, the gold nanoparticles of the plasmonic nanoantenna is divided into three layers (Figure 17): (i) one bottom layer of thickness $t_1=10$ nm interfacing the substrate, (ii) a middle layer above the first one of thickness $t_2=10$ nm, and (iii) a third top layer of thickness $t_3=20$ nm. For each case, the radiated power $P_{rad}$ is computed utilizing the corresponding volume eigencurrents [72]:

$$\begin{array}{c}\hfill P_{rad}=Re\left(P_d\right)=Re({\iiint }_{V_i}\overline{E}·\overline{J}d{V}_i)=Re\left({\iiint }_{V_i}\sigma {\overline{E}}^{2}d{V}_i\right)\end{array}$$

(25)

Explicitly, $P_d$ refers to the total dissipated power, whose real part represents the radiation power $P_{rad}$, while its imaginary part, the Joule losses $P_{loss}$.

Assuming the same input power, $P_{in}=P_{rad}+P_{loss}$, for all the three cases, the corresponding radiation efficiencies are estimated [72], ${\eta }_{R_i}=\frac{P_{ra{d}_i}}{P_{in}}$, where $i=1,2,3$ corresponds to the three different layers (bottom (red), middle (green), and top (yellow) as shown in Figure 17). Note that all dielectrics are assumed to be lossless, and any other source of radiated power other than the two gold nanoparticles is omitted.

Table 3. Radiation efficiencies and contribution to the radiated power of every layer of plasmonic nanoantenna (Figure 17) for antenna modes ($3,4$) as appear in Table 2.

	Radiation Efficiency (%)		Contribution to ${\mathit{P}}_{\mathit{rad}}$ (%)
mode	3	4	3	4
layer-1	68.26	57.22	78.64	82.86
layer-2	9.84	7.61	11.33	11.02
layer-3	8.70	4.23	10.03	6.12

lists the radiation efficiency and the corresponding contribution to the radiated power for all the different cases. One may observe that the currents flowing at the bottom layer (red) layer of the nanoparticles contribute more than $78\%$ of the radiated power, while the contribution of the remaining currents is less than $12\%$ (Table 3).

The above result is significant, as it suggests that appropriate metal thickness should be used for the design and fabrication of these antennas. However, the important point refers to the smoothness of the lower metal surface. Although the model herein considers a “perfect surface”, any roughness caused by practical fabrication may yield significant losses since the interface current will follow them. Regarding the simulation, the created mesh should ensure this perfect surface at the interface by utilizing appropriate fine discretization (i.e., 4 nm or $\lambda /160$ as shown herein). Although the thickness of metallic layers is proved non-critical, it must satisfy the necessary mechanical strength and fulfill the fabrication demands.

5. Conclusions

This work presented for the first time an in depth eigenanalysis for the study of antennas operating from terahertz to optical frequencies. Utilizing the current distributions of eigensolutions, one may investigate whether surface currents are sufficient for the correct characterization of nanoantennas, or the involvement of volume currents is needed. It is herein proved that for frequencies over 10 THz and metal thickness less than 100 nm, field penetration is prominent, requiring volume currents to be involved in the analysis. However, the examination of the radiated power contributed by the metallic layers revealed that currents penetrating inside the metallic regions have a small effect on the efficiency of nanoantennas. Hence, cost reduction can be succeeded, suggesting that a thin metallic layer, as allowed by the mechanical strength and the fabrication process, can be used. However, the high current density flowing at the metal–substrate interface asks for attention in ensuring a smooth surface during fabrication as well as fine discretization within electromagnetic simulations.