Design and Circuit Modeling of Graphene Plasmonic Nanoantennas

Nanonntennas are critical elements of nanoscale wireless communication technologies with potential to overcome some of the limitations of on-chip interconnects. In this paper, we model the optical response of graphene-based nanoantennas (GNAs) using an equivalent resistive-inductive-capacitive (RLC) circuit by incorporating plasmonic effects that are present in graphene at terahertz (THz) frequencies. The equivalent circuit model is used to estimate the resonance characteristics of the nanoantenna, thereby facilitating geometry and material optimization. Three GNA structures are considered, namely bowtie, circular, and rectangular dimers with vacuum and silicon dioxide dielectric environment. We characterize the radiation efficiency, the Purcell factor, the quality factor, and the resonant frequency of various antennas with respect to the physical properties of the surrounding dielectric media and the graphene sheet. The GNAs are designed for a resonant frequency in the range of 10 – 25 THz with field enhancements around 104 – 105 and radiation efficiency ~5 – 16%. The circuit model is applied to examine the trade-off between the radiation efficiency and the field enhancement in the antenna gap. While bowtie antennas have higher field enhancement due to charge accumulation in their pointed structure, they suffer from low radiation efficiency stemming from the large spreading resistance losses. The circuit model is validated against finite-difference time-domain (FDTD) simulations conducted in Lumerical, and excellent match between the model and numerical results is demonstrated in the THz regime. The circuit models can be readily used in a hierarchical circuit simulator to design and optimize optical nanoantennas for THz communication in high-performance computing systems.


I. INTRODUCTION
Terahertz (THz) technology has recently emerged as a highly sought-after scientific tool for a wide variety of applications [1], such as medical imaging [2], [3], biological sensing [4], [5], and molecular spectroscopy [6], in which high bandwidth performance, lightweight design, and low-power requirements are a premium. However, one of the most potentially significant applications of THz technology is in wireless communication, including networks-on-chip, [7]- [11], owing to the extremely high bandwidth of the THz band, ranging from tens to thousands of GHz depending on the signal transmission distance. Such high bandwidths of the The associate editor coordinating the review of this manuscript and approving it for publication was Mengmeng Li. THz band can alleviate the spectrum scarcity and capacity limits of current wireless technologies, which are limited to bandwidths of less than 10 GHz. An important component of the THz communication system is an optical nanoantenna (ONA) that can convert free space THz waves into localized energy and vice-versa for information processing and transmission [12], [13]. Metal-based ONAs using gold and silver show high electric field intensity enhancements [14], but the inability to dynamically tune their response over a spectral range limits their application space. The low-loss plasmonic response of graphene and its spectral tunability via electrostatic gating and geometry modification in the THz frequency range (0.1−30 THz) make it a promising candidate for ONAs in on-chip wireless communication. In particular, graphene interfacing a dielectric supports highly localized transverse magnetic (TM) surface plasmons in the THz band [15], [16]. Due to these surface plasmon oscillations, the spontaneous emission rate of a quantum emitter placed close to a graphene layer is enhanced by several orders of magnitude [17]- [20]. This enhancement in spontaneous emission, called the Purcell effect [21], can be utilized to implement efficient ONAs for on-chip light emission and detection, thus alleviating the need for traditional lasers that consume high power and area.
Several prior works have discussed the operation, design, and analysis of ONAs [12], [22]- [31]. RF engineers are equipped with the tools of circuit theory to analyze the frequency response of antennas and gain valuable insights about their resonant properties. Any resonating nanostructure can be modeled as an oscillator with radiative, inductive, and capacitive components, and the resulting radiation can be modeled as an effective radiation resistance. Circuit models for ONAs allow for application-driven design optimization and can be used to derive intuitive descriptions of the dispersion and plasmon mode propagation without invoking quantum mechanics. In prior work, the impedance of ONAs for a given spectral range has been calculated by measuring the voltage to current ratios in full-field simulations [24], [32]. A circuit representation of spherical nanostructures suspended in vacuum has been developed for operation at optical frequencies [33]. A surface integral method to derive lumped circuit elements from electromagnetic simulations has been reported for nanodipole antennas [34]. In Ref. [35], an antenna model was developed in the high-THz band (>100 THz) to study the Purcell effect in a dipole source placed in the vicinity of a spherical gold nanosphere when the system is suspended in vacuum. Circuit models have also been used to study surface plasmon propagation at THz frequencies in single and multi-layer infinite graphene sheets, nanoribbons and carbon nanotubes [33], [36]- [40]. References [41] and [42] present the design and transmission line model, respectively, for a photomixer driven radiating graphene dipole antenna (rather than simple scatterers) at very low THz frequencies. A complete treatment of modeling, fabrication, and applications of optical nanoantennas can be found in [43].
In this work, we focus on graphene-based nanoantennas (GNAs) that are characterized by their high Purcell factors, strong scattering mechanisms and sub-wavelength localization of light. A review of GNAs for THz systems can be found in [44]. The use of graphene to implement programmable metasurfaces in the THz band for potential applications in beam steering and focusing can be found in [45]- [47]. We extend prior works by developing lumped circuit models of graphene dimer structures excited by a dipole source in the THz band. The two key contributions of this work are: (1) development of an equivalent circuit model of GNAs using RF circuit theory, antenna theory and full-wave finite-difference time-domain (FDTD) simulations, and (2) characterization of GNA figures of merit with respect to the physical properties of graphene and the substrate used.
We use the classical description of the problem and atomistic details of the GNA are ignored. For a quantum-mechanical treatment of the problem based on time-dependent density functional theory and atomistic calculations, readers are referred to [48].
The remainder of the paper is structured as follows. In Section II, we present the theory of plasmonic effects in graphene and their use toward generating and detecting THz radiation in various GNA configurations. The equivalent circuit model of the GNA is presented in Section III, while physical equations of various circuit components are derived in Section IV. The effect of changing material parameters of graphene on the resonant features of various GNA structures is also analyzed in Section IV. The equivalent circuit model of the GNA is validated against full-field electromagnetic simulations in Section V. This section also details the differences in the resonant characteristics of various dimer structures. Finally, the paper concludes in Section VI.

II. OPERATION OF THE GRAPHENE NANOANTENNA (GNA)
GNAs are an emerging field of optics and utilize plasmonic effects in graphene to localize optical energy below the diffraction limit of light [49]. Unlike conventional RF antennas for which the resonant length is equal to half the wavelength of incident light, GNAs can be designed to have a much shorter resonant length [50]. Accurate modeling of GNAs at THz frequencies requires a detailed understanding of polarization and displacement currents and the finite skin depth associated with surface current conduction [13]. In GNA plasmonic structures, a constructive interference of multiple scattered or re-radiated fields can cause local field enhancement; therefore, nanostructures that can maximize the total effective area that the incident field perceives when interacting with the nanoantenna are highly desired. The interaction of the incident radiation with the nanoantenna is quantified through the extinction cross-section (α ext ), defined as the sum of the scattering cross-section (α scat ) and the absorption cross-section (α abs ). The scattering cross-section gives the ratio of the scattered/re-radiated power to the incident power flow per unit area, while the absorption cross-section gives the ratio of absorbed power to the incident power flow per unit area. In Ref. [29], it was shown that the interaction of THz radiation with a graphene nano-patch antenna with ∼10's of µm 2 cross-sectional area is dominated by absorption.
We can translate the classical concept of antenna aperture, which measures the efficiency of interaction between the incident radiation and the antenna, to the THz regime by defining an ''interaction'' cross-section of the dipole source coupled to the nanoantenna. In the next sub-section, we discuss the localization of optical energy in graphene nanostructures due to plasmonic effects. Overall, our results shed light on the optimal design space of GNAs and can be used as guidelines to augment the fabrication and characterization of GNAs.

A. SURFACE PLASMONS IN GRAPHENE AND LOCAL FIELD ENHANCEMENT
It is well-established that graphene placed in a dielectric media can support TM surface plasmon modes at THz frequencies [51]. The TM plasmon dispersion relationship, k sp (ω), in graphene for plasmons propagating along x− direction as shown in Fig. 2 is found from where i = √ −1 is the imaginary number, ω is the frequency of operation, c is the speed of light in vacuum, 0 = 8.85 × 10 −12 F/m is the permittivity of free space, and 1 and 2 are the relative dielectric permittivities of the media surrounding the graphene sheet. Considering only intra-band scatterings in the THz regime, the dynamical conductivity of graphene, σ (ω), using the Kubo formalism is given as [52], [53] Here, e is the elementary charge,h is the reduced Planck's constant, T is the temperature, µ cp is the electrochemical potential of the graphene sheet, and f d = is the Fermi-Dirac statistics (k is Boltzmann constant). The factor τ is a phenomenological energy-independent electron relaxation time ( = τ −1 is the electron relaxation rate) that accounts for losses due to electron-impurity, electron-defect, and electron-phonon scatterings. Carrying out the integral in (2) and assuming N layer number of decoupled parallel layers in the multi-layer graphene stack, we obtain Since N layer and µ cp appear as multiplicative factors, one could adjust the value of χ by changing either N layer or µ cp . In this work, we choose χ = 1.4×10 11 ( s) −1 , which can be obtained for N layer = 1 (N layer = 5) and µ cp = 1.2 eV (µ cp = 0.24 eV). Additionally, the effect of electrostatic screening in the multi-layer graphene structure can be included by adjusting the effective number of layers in (3b). Note that model in (3) ignores spatial dispersion effects and is strictly applicable to multi-layer graphene fabricated via chemical vapor deposition technique in which the layers are decoupled and the linear dispersion relationship with a Dirac cone is preserved [54], [55]. Since, ω < ω inter (inter-band threshold frequency) and ω < ω op = 0.2 eV (optical phonon threshold), plasmon decay in graphene is  to be dominated by elastic scatterings via acoustic phonons (AC) and charged impurities (CI) [56]. Using Matthiessen's rule, the total scattering rate in graphene is given as [57] = AC + CI , where D ac is the acoustic deformation potential, ρ m is the mass density, v ph is the phonon velocity in graphene, Z is the net charge of the impurity, avg = 1 2 ( 1 + 2 ) is the average permittivity of the dielectric media surrounding graphene, N CI is the density of charged impurities, and N s is the carrier density in graphene. The value of N s is obtained from the energy-dependent integral of the density of states, including the disorder due to electron-hole puddles around the Dirac point in graphene, and the Fermi-Dirac statistics as explained in [58]. Figure 1 shows the room-temperature scattering rate versus µ cp in graphene for different values of N CI , and the inset shows the mean free path of electrons in graphene versus µ cp . While the scattering rate for phonon-dominated collisions increases with an increase in µ cp , for collisions due to charged impurities, scattering rate decreases with an increase in µ cp . This behavior is also observed from (4). In a pure graphene sample, is limited due to the acoustic phonon scatterings. For D ac = 7.25 eV, µ cp = 1.2 eV, and N CI = 2 × 10 14 cm −2 , we get = 8.5 × 10 11 s −1 , which is the value chosen for all simulations reported in this paper, unless otherwise noted. The mean free path of electrons corresponding to = 8.5 × 10 11 s −1 is 1.17 µm at 300 K. Similar values of electron mean free path have been reported experimentally in DC transport measurements in graphene [59]- [61]. Other scattering mechanisms, such as due to optical phonon emission, surface polar phonons of the substrate, and short-range scatterers (adatoms adsorbed on graphene, vacancies) if present can be included in as discussed in [57]. Nevertheless, could essentially be treated as a phenomenological parameter that is extracted based on experimentally measured mobility in the sample.
The material parameters of graphene, and χ, are fixed with respect to frequency and the dielectric environment in FDTD simulations and are used for model validation in Sec. V. We note that the specific choice of and χ values does not affect the development of the equivalent circuit model of the GNA and the physical insight gained from it. A different value of material parameters would, however, yield different numerical values of the GNA circuit components and its figures of merit.
We use (3) in (1) in the non-retarded regime (i.e. q ω/c) to re-write the plasmon wave-vector as When the frequency of the incident radiation matches the frequency of the induced dipole moment in the sample, a peak in α scat and α abs is observed. The maximum value of the peak and the plasmon linewidth are limited by ohmic losses (scattering rate) in the material and the radiation damping [62]. The wavelength of the plasmons (λ sp ) is smaller than the wavelength of the incident radiation (λ light ). For example, at a frequency of 25 THz, the ratio of plasmon wavelength in a graphene sheet suspended in vacuum to the wavelength of the incident light is ≈0.17 for χ = 1.4 × 10 11 ( s) −1 . This implies that the energy contained in the incident wave is squeezed into a smaller volume resulting in considerable local enhancement of the electric field intensity. The plasmon confinement is higher for graphene sheets with a lower electrochemical potential and for substrates with a high dielectric permittivity. Moreover, the degree of electric field enhancement due to plasmonic effects depends on the details of the charge distribution in graphene. Two cases for a graphene sample placed under radiation are illustrated in Fig. 2. A truncated graphene sheet, shown in Fig. 2a, can radiate in the THz band due to its broken translational symmetry. However, in this case, it is the fundamental plasmon mode with sinusoidal charge distribution in the plane of the graphene sheet that couples to the incident radiation [63]. Dimer structures, on the other hand, allow higher order modes to couple together and cause very high electric field enhancements in the gap. As illustrated in Fig. 2b, the charge distribution in dimer structures is intensified at the gap. As a result, the spontaneous emission rate of a quantum emitter placed in the gap of a dimer structure can be enhanced, making it a suitable candidate for implementing a GNA.   The various GNA structures analyzed in this work are illustrated in Fig. 3. In all structures, the GNA comprises two single layer graphene arms of equal dimensions separated by a gap of distance g and placed on a substrate. The source of excitation is a dipole with a length of 3 nm placed in the gap of the GNA. The substrates considered are vacuum and silicon dioxide (SiO 2 ), which is a commonly used substrate for graphene deposition.

III. EQUIVALENT CIRCUIT MODEL
In this section, we present an equivalent RLC circuit model for the GNA structures shown in Fig. 3, while the antenna dimensions are shown in Fig. 4. The equivalent circuit model VOLUME 8, 2020 is used to characterize the figures of merit of the GNAs, namely the resonant frequency, the radiation efficiency, and the Purcell factor. The Purcell effect refers to the modification of a quantum system's spontaneous emission rate by its electromagnetic environment. Despite the quantum mechanical origin of the Purcell effect, the underlying electrodynamics of an GNA is similar to classical antennas due to the linearity of Maxwell's equations, thus allowing their emission and resonant properties to be modeled using the classical circuit theory. The gap distances considered in this paper are large enough (>5 nm) to ignore the effects of quantum tunneling between the antenna arms, thus alleviating the need for a quantum mechanical treatment [24]. Finally, we note that the the circuit model derived here is applicable over a broad range of THz spectrum, up to frequencies below the inter-band threshold. In this paper, we specifically design the antennas such that their resonant frequency is in the range of (10-25) THz for all configurations. By appropriately choosing the dimensions, dimer configuration, and electrostatic gating, a GNA can be optimized to radiate at different frequencies within the THz band, as discussed briefly in Sec. IVB. Readers may refer to previously published works [29], [41], [50] for a deeper analysis of GNA design and performance at frequencies less than 10 THz.

A. CIRCUIT MODEL WITHOUT THE GNA
The circuit model of a radiating dipole source without a GNA is shown in Fig. 5. The quantum mechanical current associated with a dipole source radiating at a frequency of ω is given as where I 0 = q 0 ω is the peak quantum mechanical current. For a small dipole length, the current can be considered to have a constant phase along the length of the dipole and the quantum mechanical current can be approximated as I = qω. The total decay rate of this radiating dipole source, also called a quantum emitter, in the absence of the GNA is the sum of the radiative (γ r0 ) and the non-radiative (γ nr0 ) decay rates. That is, γ 0 = γ r0 + γ nr0 . In terms of the power radiated, γ 0 = P 0 /hω, where the total power is given by P 0 = P r0 + P nr0 [64]. The efficiency of the emitter in a homogeneous environment is given as In terms of the resistance of the quantum emitter, the efficiency can also be expressed as where R r0 is the radiation resistance of the dipole and R nr0 is the non-radiative loss resistance of the dipole source in the absence of the GNA. The radiation resistance R r0 = π/6(2x d /λ) 2 Z 0 n, where x d is the dipole length, λ is the incident radiation wavelength, Z 0 (= 377 ) is the impedance Here, R r0 and R nr0 are the radiative and non-radiative resistances, respectively.
of free space, and n is the average refractive index of the surrounding media. Intrinsic non-radiation loss mechanisms vary for different emitters and include, for example, phononscattering, Auger recombination of multiexcitons in colloidal quantum dots, and losses due to lattice defects in colored centers. For a dipole with unity quantum yield, Q 0 = 1, the intrinsic non-radiative decay is ignored.

B. CIRCUIT MODEL IN THE PRESENCE OF THE GNA
The equivalent circuit diagram for the GNA with a quantum emitter in the gap of its arms is shown in Fig. 6. The capacitive elements in this circuit model include C ant , which is the capacitance due to the charges on the antenna arms, while C couple denotes the coupling capacitance of the arms. The coupling capacitance is due to the sidewall and fringing electric field lines that exist between the antenna arms. Since the thickness of graphene is insignificant compared to the gap distance for all GNAs considered here, most of the coupling capacitance is due to the fringing capacitance (C top−bot couple ), the origin of which is the electric field lines extending between the top and bottom surfaces of the GNA arms. While the sidewall coupling capacitance (C side couple ) scales with the ratio t/g, where t is the thickness and g is the dimer gap of the GNA, the top and bottom surface capacitance scale with the ratio L/g, where L is the length of the GNA arms. The substrate height, chosen as 2 µm for all designs, is significantly greater than the GNA gap distance (≈25-30 nm for all structures), which leads to negligible capacitance associated with the ground-plane electric flux. There are two inductive components in the model: L K is the kinetic inductance due to the inertia of electrons (plasmonic effects), and L f is the Faraday inductance. The ohmic losses in the antenna are modeled using the resistance R , while the radiation is represented using R r . The radiation resistance includes both the dipole source radiation resistance and the extra radiation resistance due to the presence of the GNA that leads to enhanced emission. That is, R r > R r0 . The ohmic resistance of the graphene sheet depends on its shape, carrier density, and the physics of carrier scattering.

1) RESONANT FREQUENCY
For given LC parameters, the resonant frequency of the GNA is given as 129566 VOLUME 8, 2020 FIGURE 6. Circuit diagram of a dipole current source, I QM , with a nanoantenna. The resistance includes radiation resistance, R r , and the non-radiative ohmic resistance, R . The capacitance comprises the coupling capacitance, C couple , and the antenna capacitance, C ant . The inductive components of the antenna are due to the Faraday inductance, L f , and the kinetic inductance, L K . The current flowing through antenna arms is I A , which is enhanced compared to I QM for a greater than unity quality factor and when coupling capacitance is much smaller than antenna capacitance.
where C couple = C side couple + C top−bot couple . In terms of the ratios C top−bot couple /C ant = α top−bot and C side couple /C ant = α side , the net capacitance is simplified to The dependence of the resonant frequency, ω res , on the ratio α top−bot is shown in Fig. 7. For sidewall capacitances that can be neglected in comparison to the antenna capacitance, the resonant frequency decreases as α top−bot increases and it saturates for α top−bot 1. In the large-α top−bot limit, ω res is limited by the antenna capacitance (C net ≈ C ant ). In this limit, the large coupling capacitance shorts the gap of the dimer GNA structure such that the resonant frequency converges to that of a simple wire or patch GNA. Likewise, a large value of α side lowers the resonant frequency making the dimer structure equivalent to a wire or patch antenna.

2) MAXIMUM RADIATED POWER
The antenna arms modify the quantum mechanical current of the emitter, which is given as [19] where x d is the dipole length and g is the gap length of the GNA. A simple circuit analysis based on Fig. 6 shows that the current flowing through the GNA arms at resonance condition is given as where Q f is the quality factor of the GNA and is a measure of the antenna bandwidth (BW). That is, Q f = ω res /BW. Using I A from (12), the maximum power radiated by the GNA at resonance is given as For

3) ANTENNA EFFICIENCY
The antenna efficiency depends on resistive elements of the circuit and is given as where γ r and γ nr are the radiative and non-radiative decay rates, respectively, of the emitter in the presence of the GNA, while γ is the ohmic loss decay rate in the GNA. We use (7) in (14) and assuming γ nr = γ nr0 , we obtain: In the THz band, λ x d (x d = 3 nm is the dipole length), implying R r0 ∼ 10 −4 . Given that the radiation and loss resistances of the GNA are orders of magnitude larger than R r0 , it is reasonable to neglect the term R r0 (1−Q 0 )/Q 0 in (15), unless Q 0 10 −3 . Here, we model the radiation efficiency of the GNA using the following: In addition to ohmic loss in the antenna, there will be a dipole-associated spreading resistance loss that is not inherent to the fundamental GNA mode. The spreading resistance, R spread , appears in series with the current source I QM in Fig. 6.

VOLUME 8, 2020
The antenna efficiency in the presence of spreading resistance loss is modeled as The above equation reveals that the effect of the spreading resistance in determining the antenna efficiency is suppressed by a factor of (I QM /I A ) 2 ≈ 1/Q 2 f .

4) PURCELL FACTOR
The Purcell factor is defined as the ratio of the power emitted in the antenna (including antenna losses) to the power emitted by the dipole in a homogeneous environment (without the antenna). That is, In terms of the antenna radiation efficiency, Q ant , and the quantum efficiency of the emitter, Q 0 , which is further simplified to for unity quantum efficiency of the emitter. Here, P r and P refer to the power radiated and ohmic power loss, respectively, in the antenna. The enhancement of spontaneous emission (E) of the quantum emitter due to the GNA is defined as P r /P r0 , where P r0 = (qω) 2 R r0 is the radiative power of the emitter. Using P r from (13), we get where c is the speed of light in vacuum. The enhancement factor, E, can be increased by reducing the gap between the GNA arms and reducing the non-radiative losses in the GNA. Once the enhancement factor is known, the Purcell factor of the GNA is evaluated as

IV. MODELING OF CIRCUIT ELEMENTS A. CAPACITIVE ELEMENTS
The sidewall coupling capacitance of the GNA is given as where A tip is the cross-sectional area of the GNA tip. For the bow-tie GNA, A tip = w t, where t = (2N layer − 1) t (here t = 0.35 nm is the thickness of each graphene layer and N layer = 10), and w is shown in Fig. 4. In the case of rectangular GNA, A tip = Wt, while for the circular GNA, A tip = πrt/2. Considering the small thickness of the graphene sheet and the small tip area, C side couple C ant for all GNAs considered in this work. The top-and bottom-wall coupling capacitance, C top−bot couple = α top−bot C ant , scale with the width of the GNA and depend on the ratio of the GNA length and the gap distance. Since the GNA length is three orders of magnitude larger than the GNA thickness, we expect C couple ≈ C top−bot couple . The value of the coupling capacitance is extracted from FDTD simulations for each antenna geometry.
The antenna capacitance is modeled as where x− (y−) direction is defined along the length (width) of the GNA, and n s (x) is the two-dimensional (2D) charge distribution in the antenna arms. For uniform charge distribution in the 2D plane of the GNA, C ant evaluates to where A GNA is the area of the GNA. However, in dimer configurations, the symmetry of the structure is broken and the charge distribution is non-uniform as illustrated in Fig. 2. This charge distribution is of the form where the parameters n 0 and x 0 are empirical parameters that depend on the GNA configuration and dimensions. Here, we incorporate the effect of non-uniform n s (x) by multiplying C ant with an empirical factor k c , which accounts for the fact that the surface charge distribution increases rapidly from the center of the arm toward the antenna gap. The charges that accumulate in the corners of the antenna arms add to the capacitance, implying k c > 1. This factor is extracted from numerical FDTD simulations and is a constant, in the range of (1.5 − 2.5), for a particular configuration. While close to the Dirac point, intrinsic quantum capacitance of the graphene sheet must be included in the circuit model, here we ignore it since µ cp (= 1.2 eV) φ t for all simulations.

B. INDUCTIVE ELEMENTS
The oscillating surface charges in the nanoantenna constitute electric current and, therefore, produce a magnetic field around the graphene sheet. This effect can be captured by defining a Faraday inductance associated with the nanoantenna, which is given as where L is as shown in Fig. 4, k sp is defined in (1), and the factor ζ = 2r for the circular nanoantenna, and ζ = W for the rectangular and bow-tie nanoantennas. In the non-retarded regime (k sp k 0 ) assumed in this paper, the plasmon wavelength is quite short. A consequence of this is that the Faraday inductance becomes insignificant in comparison to the kinetic inductance [37], and, therefore, has no impact on the resonant features of the GNA. Indeed, we confirm numerically and analytically that the Faraday inductance of all GNAs studied here is less than 3.5% of their kinetic inductance.
The kinetic inductance of the GNA is given as where the function g(x) depends on the GNA geometry. For rectangular GNA, g(x) = W , and for circular GNA, g(x) = 2r. However, due to non-uniform width of the GNA arms in the bowtie geometry, g(x) is given as The integration of g(x) in (28) is Based on the dimensions shown for various GNAs in Fig. 4, G c = 3.88 and 1.75 for bowtie and rectangular GNAs, respectively. In (28), σ r (σ i ) is the real (imaginary) component of the 2D sheet conductivity of graphene obtained at the resonant frequency. Using Eqs. (3a) and (3b) Per the above equation, the value of L K is dependent on the dielectric media surrounding graphene through the factor /ω res . For = 0.85 THz and f res ≥ 10 THz, L K ≈ G c /χ and is independent of the dielectric media surrounding graphene (given fixed χ and ).

C. RESISTIVE ELEMENTS
The loss resistance of the GNA comprises the ohmic resistance and the spreading resistance. The total loss resistance of graphene is given as Since and χ are chosen as fixed parameters, the difference in R across various GNA configurations originates from the geometric correction G c . The spreading resistance, R spread ∝ 1/σ r (ω), where the constant of proportionality is a strong function of the GNA configuration and the dielectric environment. To estimate the value of R spread , we could treat the induced currents on the antenna tips as a point contact of size βt, where t is the thickness of the graphene sheet, and β is a factor that is related to the geometry of the antenna tips [64]. Our FDTD simulations show that β is in the range of (0.03-0.05) for various GNA structures and is influenced strongly by the distribution of induced current in the antenna arms. We also find that R spread ∼ (20 − 35) × R and is higher for bowtie and rectangular GNAs due to the higher confinement of electric field lines in the GNA arms around the gap. As such, the total loss resistance for GNAs is dominated by their spreading resistance. However, the impact of R spread on Q ant is suppressed by a factor of Q f , which in our geometries is much greater than unity due to the narrowband emission characteristics of the GNAs. With increased dielectric permittivity of the media surrounding the GNA, the impact of spreading resistance with respect to ohmic resistance is reduced for all GNA types. Due to the plasmonic effects in graphene, the effective wavelength λ eff is much shorter than the incident wavelength, λ, and is empirically given as λ eff = n 1 + n 2 λ/λ sp , where n 1 and n 2 are geometry-and environment-dependent parameters [31]. For a thin-wire antenna at RF frequencies where plasmonic effects are negligible, the radiation resistance is given as R r = 30π 2 (L/λ) 2 . However, the presence of plasmonic effects can be accounted for in the radiation resistance by replacing λ with λ eff . Using the concept of effective wavelength, the radiation resistance of half-wavelength GNA is given as [13] Since λ eff λ, the radiation resistance of GNAs with plasmonic effects is only around few 's. For GNAs analyzed in this work, λ eff /λ is ≈(0.2-0.3), while λ eff /λ sp is in the range of 1-2.

D. IMPACT OF MATERIAL PARAMETERS
The circuit components in Fig. 6 depend on the electrochemical potential in graphene as well as the shape and size of the dimer structure. Here, we demonstrate the impact of adjusting µ cp on the antenna's resonant features while keeping its size and, therefore, its geometric correction factors the same. The size of various GNA structures is given in Fig. 4.
The resonant frequency of the GNA is given as ω res = 2L k α top−bot υ −1/4 , where υ = 2 0 avg 2 A GNA k c . As avg increases, ω res decreases as −1/2 avg ; this scaling trend is corroborated by our FDTD simulation results presented in Sec. V. Figure 8 shows the GNA resonant frequency versus µ cp for all dimer structures. Since L K and C net both decrease linearly with µ cp , ω res varies as √ µ cp for all GNA configurations. For µ cp = 1.2 eV and the dimensions chosen in Fig. 4, the resonant frequency is greater than 20 THz for all GNAs. Note that does not have an impact on the resonant frequency unless it is comparable to ω res as shown in the inset of Fig. 8.
The effect of electrochemical potential on antenna efficiency is analyzed in Fig. 9(a) for various GNA structures.
In the presence of a finite spreading resistance, the antenna efficiency shows a non-monotonic behavior with respect to µ cp : a maximum in Q ant occurs at µ cp ≈ 0.75 eV. If the spreading resistance were absent, the antenna efficiency would increase by a factor of 3× for rectangular and bowtie GNAs, while the improvement in the efficiency of circular GNAs would be around 1.25×. This is an expected because spreading resistance is quite significant in pointed GNA structures as discussed earlier and validated using FDTD simulations in Sec. V. Figure 9(b) shows the power radiated by the GNA at resonance as a function of its electrochemical potential. The increase in P rad versus µ cp is quite rapid when µ cp is below 0.5 eV. However, P rad shows quasi-saturation with µ cp because eventually the loss resistances in the GNA dominate its net resistance, limiting the power the antenna can radiate.

A. MODEL VALIDATION AGAINST FDTD
FDTD simulations are performed for each GNA configuration to study its response in the THz band. The FDTD, in principle, is a full-wave electromagnetic modeling scheme that accurately captures optical modes in photonic/optical systems through a numerical solution of the Maxwell's equations. Here, we use the Lumerical package to carry out the FDTD routine. The numerical simulation utilizes the differential form of Maxwell's equations with a time-stepping algorithm that approximates all derivatives in Maxwell's equation. The spatial grid in FDTD is a Yee-lattice in which electric and magnetic fields are computed on separate grids and at successive half-steps in time. The FDTD is a time-domain method and therefore applicable for generating broadband antenna response. In Fig. 10, we show the total, far-field, and absorbed power of the GNA suspended in vacuum and deposited on SiO 2 . The power is defined as W (ω) = S(ω) · d , where defines the surface area over which the Poynting vector S(ω) is integrated. The Poynting vector is constructed from the Fourier transforms of the time-dependent electric fields [65]. The GNA power clearly shows the presence of side lobes, and this power-frequency profile can be described analytically as P r = P 0 sinc 2 ((ω − ω res )/BW), where the GNA bandwidth, BW, depends on the quality factor, and P 0 is the maximum power radiated at the resonant frequency. From these simulations, we measure the following figures of merit for each antenna configuration: resonant frequency, radiation efficiency at the resonant frequency, and the Purcell factor. Results are noted in Table 1. For a fixed ohmic loss in the GNA sheet, the response of the GNA degrades as the relative dielectric permittivity of the substrate increases. Yet, for all cases examined here, the resonant frequency falls within the expected range of 10 THz to 25 THz. The Purcell factor of each GNA configuration stays well above 10 3 with a radiation efficiency in the range of ∼(5-16)%. One viable route toward improving the radiation efficiency for all GNA configurations is to reduce the sources of scatterings in graphene by optimizing the deposition process or by using dielectrics that are lattice matched to graphene.
Next, we validate the figure of merit of the GNA obtained from the equivalent RLC model against numerical results reported in Table 1. The figures of merit evaluated in our work  include the maximum radiated power, antenna quality factor, the Purcell factor, and the maximum antenna efficiency when the resonant condition is attained in the antenna, i.e. ω = ω res . By combining analytic models with FDTD simulations, we obtain the circuit parameters of the GNA comprising seven components: R r , R , R spread L K , L f , C ant , C top−bot couple . Since our circuit model has seven elements, one can utilize seven measurements from an FDTD simulation to obtain their values. While such an approach has been adopted in previous works [19], [34], [64], it lacks in its capability to predict antenna response by tuning the geometry and material parameters of the antenna. On the other hand, our goal in this work is to link the physics and design of the GNA to its performance characteristics in the THz band. Therefore, out of the seven unknown circuit elements, we model four of them, namely R , L K , L f , C ant , based on the material specifications of the GNA. We only use three measurements from an FDTD result, namely the resonant frequency, the resonant power and the quality factor, to extract the value of R r , R sp , and α top−bot for all GNA structures. The extracted value of R r matches very well with the physics described in (32) and further establishes a connection between circuit parameters and the physics of GNA resonance. Our rationale to use FDTD results is to arrive at an optimal trade-off between complexity and accuracy of circuit modeling. Hence, only a limited number of parameters in the model are extracted from simulation data, and not directly modeled from first principles. By combining FDTD simulations with analytic equations of various model parameters, our approach overcomes the limits of previous works where all circuit elements were retrieved from electromagnetic simulations.
The radiation resistance is extracted as while the spreading resistance is extracted as   . The factor k c = 1.5 (circle), 2.5 (bowtie), and 1.7 (rectangle). The absolute percentage deviation between a modeled (y mo ) and its measured (y me ) quantity is calculated as |y mo − y me |/y me × 100. There is negligible error in max. P r , Q f , and ω res , as these FDTD measurements are utilized to extract the value of coupling capacitance of the antenna and its spreading resistance.
Here the superscripts {ex, me, and mo} stand for extracted, measurement, and modeled values, respectively. While the emitter quantum efficiency, Q 0 , does not impact Q ant in the THz band, it is nonetheless important for the Purcell factor. From the measurement of F p and Q ant via FDTD simulations, the efficiency of the quantum emitter is extracted as Q ex 0 = F me p Q me ant /E mo . In our approach, we need to conduct FDTD simulations only once at the resonant frequency to obtain the model parameters, while the GNA response at any other frequency is obtained using the circuit model of Fig. 6.
Model parameters for each GNA are listed in Table 2. The figures of merit of the GNA obtained using the equivalent circuit model are compared against those obtained from numerical simulation as shown in Figs. 11(a) and (b). The model does an excellent job at capturing the key metrics of the GNA at the resonant frequency with an error of less than 35% in Q ant and F p and negligible error in Q f , ω res , and peak P r in all cases. Given the simplicity of the circuit model and the intuitive insight it provides, the error is reasonable. Using FDTD simulation results to extract the numerical values of all circuit elements would lead to a lower error rate, although this approach would not be able to capture the impact of materials on the GNA performance.
Ignoring the spreading resistance overestimates the value of Q ant , particularly for bowtie and rectangular GNA configurations. This is due to the high field confinement in these structures which increases the contribution of spreading resistance to the net loss resistance. While the value of ohmic resistance is independent of the dielectric permittivity as evident from the relationship in (31), both radiation resistance and spreading loss resistance depend strongly on the physical properties of the dielectric. The reduction in radiation resistance with an increase in dielectric permittivity of the substrate leads to a decrease in the antenna efficiency for high permittivity substrates. In Figs. 11(c) and (d), the response of the GNA obtained from equivalent circuit model is compared against FDTD results over a broad frequency spectrum. The model faithfully captures the key features of the GNA response and accurately predicts the maximum output power at resonance, quality factor, as well as the resonant frequency. As a result of its implementation to trace the main lobe response, the model cannot capture the side lobes accurately.
The model estimates for Purcell factors, assuming 100% efficiency of the quantum emitter, are (5 − 25)× greater than the results obtained from FDTD simulations. However, the efficiency of the quantum emitter Q 0 20% for all GNAs examined here. For these values of Q 0 and R r0 < 10 −4 , the term R r0 (1 − Q 0 )/Q 0 is at most a few milli-ohms. Therefore, the approximation in arriving at (16) from (15) is justified.

VI. CONCLUSIONS
Graphene nanoantennas (GNAs) operating in the THz band can enable high-throughput wireless links at the nanoscale. Such wireless links can potentially overcome the losses associated with on-chip waveguide interconnects. To facilitate the design and optimization of on-chip wireless links, we develop an equivalent circuit model that captures the THz response of GNAs using simple resistive-inductive-capacitive (RLC) elements. Three geometries of the GNA, namely circular, bowtie, and rectangular with a gap excited by a quantum emitter, are modeled. The circuit model comprehends ohmic and spreading resistance loss as well as the effect of kinetic inductance originating from the presence of transverse magnetic surface plasmons in the graphene sheet. Based on the physics of plasmon confinement and propagation, the circuit models are derived in terms of the 2D carrier density and terahertz conductivity of the graphene sheet. Moreover, the effect of geometry and electric field enhancement at the tip of specific GNAs is also included. The circuit model is validated VOLUME 8, 2020 by comparing the results of resonant frequency and antenna efficiency against FDTD simulation results. The model provides physical insight into the operation of the GNA and can be used readily by circuit designers for technology-circuit co-design and optimization. Future work will explore power transmission of the GNA in the far-field and comparing its performance against that of plasmonic waveguides such as metal-insulator-metal and strip lines.