Plasmonic excitations in graphene nanostructures provide a particularly effective means to enhance light–matter interactions at THz frequencies. Here, we investigate the use of graphene nanoribbons for narrowband THz light emission based on the excitation of plasmonic oscillations under current injection and their resonant decay into free-space radiation. A detailed theoretical model of the underlying plasmon-enhanced thermal emission mechanism is presented, whose predictions are in good agreement with the recent experimental demonstration of this phenomenon. This model highlights the key role played by the nanostructure absorption efficiency to maximize the output radiation at the plasmonic resonance frequency. Based on this idea, we explore the integration of graphene nanoribbons with nearby metallic antennas in an open cavity configuration in order to promote critical coupling to free-space radiation and correspondingly enhance the absorption (and, therefore, radiation) efficiency by up to two orders of magnitude. The simulation results indicate that this approach is promising for the development of novel THz sources with technologically relevant emission characteristics.

## I. INTRODUCTION

Plasmonic excitations in graphene films and nanostructures offer several desirable features for fundamental studies and device applications, including a particularly strong sub-wavelength confinement, relatively large propagation lengths, and dynamic tunability through the application of a gate voltage.^{1–5} In addition, these resonances typically occur at mid-infrared and terahertz wavelengths that are not immediately accessible with more traditional plasmonic systems based on noble metals. The ability to cover the THz spectral region is particularly interesting from a technological standpoint,^{6} because THz optoelectronics still suffer from the lack of practical devices that can meet existing application needs (e.g., for THz spectroscopy, imaging, and sensing). Following extensive studies of their basic physical properties,^{7–14} initial device applications of graphene plasmon polaritons (GPPs) at THz frequencies have already been demonstrated, including photodetectors^{15–18} and optical modulators.^{19,20} Plasmon-enhanced interband light emission under conditions of population inversion has also been proposed for the development of graphene THz lasers.^{21–23}

More recently, narrowband THz light emission from graphene nanoribbons has also been reported based on the generation of hot carriers under current injection and their subsequent energy relaxation through the excitation of GPPs.^{24} Free-space THz radiation is then emitted by the resulting collective oscillations of the graphene electron gas in the nanoribbons. A similar radiation mechanism was originally investigated in early work with traditional semiconductor heterojunctions such as Si/SiO_{2} and GaAs/AlGaAs.^{25–27} Pronounced THz emission peaks were measured with these systems coupled to a diffraction grating, but only at cryogenic temperatures. Graphene can generally be expected to be particularly well suited for this device application because of several distinctive properties. First, by virtue of its low electronic heat capacity and weak electron–phonon coupling,^{28} graphene is an ideal material system for the excitation of hot carriers. In fact, prior work has already established its promise as a stable and efficient thermal emitter at mid-infrared^{29,30} and visible^{31} wavelengths. Second, the particularly strong confinement of GPPs in the supporting electron gas is favorable to enhance all sorts of light–matter interactions, including plasmon generation via hot-carrier decay.^{32} Finally, the relatively low damping rates of GPPs are advantageous to produce proportionally narrow plasmonic emission spectra even under ambient conditions, especially in high-mobility samples.^{33}

Nevertheless, in the initial report of Ref. 24, plasmonic THz light emission from graphene could only be detected up to a maximum substrate temperature T_{base} of about 200 K, with very low wall-plug efficiency on the order of a few 10^{−8}. In this work, we present a detailed theoretical model of the underlying radiation process, which can be described as thermal emission resonantly enhanced by the excitation of graphene plasmons. The predictions of this model are found to be in good agreement with the measurement results of Ref. 24 and are then used to identify the key limiting factors of the initial experimental devices. This analysis highlights several possible avenues toward improved emitter performance, including engineering the dielectric environment of the graphene plasmonic oscillators in order to promote critical coupling to external THz radiation. Following these predictions, we develop a specific design involving graphene nanoribbons simultaneously coupled to an optical cavity and THz antennas. Our simulation results show that this design can enable room-temperature operation with a two-order-of-magnitude increase in wall-plug efficiency compared to the devices of Ref. 24. Correspondingly, technologically significant THz output powers approaching the mW range become accessible with large but realistic input electrical powers on the order of 100 W delivered to a total graphene area of a few mm^{2}.

## II. PLASMON-ENHANCED GRAPHENE THERMAL EMISSION MODEL

The basic device geometry used for the initial demonstration of plasmonic THz light emission from graphene^{24} consists of multiple graphene nanoribbons connected between the source and drain contacts in a back-gated field-effect-transistor (FET) configuration [Fig. 1(a)]. The graphene active material was grown by chemical vapor deposition (CVD) on a Cu film and then transferred on the Si/SiO_{2} FET substrate with a poly(methyl methacrylate) (PMMA) support polymer. For practical convenience, perpendicular “bridge” sections [not shown in Fig. 1(a)] were also introduced in the patterned graphene samples to connect neighboring nanoribbons to one another, so as to minimize the impact of material defects and cracks on the current flow through the entire device. Graphene nanoribbons in general support localized plasmonic resonances associated with plasma waves bouncing back and forth between the ribbon edges. The resulting resonance frequencies ν_{pl} are determined by the phase condition requiring constructive interference after each roundtrip and can be tuned by varying the nanoribbon width *w* and the graphene charge density N according to a simple approximate analytical formula, ν_{pl} ∝ N^{1/4}/*w*^{1/2}. These localized plasmonic excitations have been studied extensively in recent years in several transmission spectroscopy measurements at both mid-infrared^{34,35} and THz^{7,10} wavelengths.

A representative THz emission spectrum measured with this graphene nanoribbon device geometry is shown in Fig. 1(b).^{24} The specific sample used in these measurements consists of 530-nm-wide nanoribbons arranged in a 150 × 200-*μ*m^{2} array with 1-*μ*m pitch, held at a substrate temperature T_{base} of 80 K with a carrier density N of 5.4 × 10^{12} cm^{−2}, and driven with an input electrical power P_{in} = I_{DS}V_{DS} of 0.2 W (here I_{DS} and V_{DS} are the applied source-drain current and voltage, respectively). The trace plotted in Fig. 1(b) is normalized to the emission spectrum of the same sample gated at charge neutrality under otherwise identical conditions, so as to fully highlight the plasmonic contribution. A sharp peak with a quality factor of about 4 is observed at a center frequency of 7.9 THz, in excellent agreement with predictions based on the above-mentioned analytical expression for ν_{pl} as well as numerical simulations based on the finite difference time domain (FDTD) method.^{24} The output power density per unit frequency per unit graphene-sample area Δρ_{out} quantified on the vertical axis of Fig. 1(b) was estimated by comparing the measured strength of the plasmonic radiation peak with the graybody emission background of the same device.^{24} Based on this estimate, we obtain a wall-plug efficiency ΔP_{out}/P_{in} of 1.7 × 10^{−8}, where ΔP_{out} is the total plasmonic output power integrated over all frequencies and across the entire device area.

The radiation process of Fig. 1(b) involves current injection in the nanoribbons to heat up the graphene electron gas, leading to an increased population of GPPs and therefore enhanced thermal emission at the nanoribbon plasmonic resonance frequency. This process can be described using a general theoretical model for thermal emission from a collection of resonant nanoparticles.^{36} In this model, consistent with Kirchhoff's law of thermal radiation, the spectral radiance S_{out}(ν, T) (output power per unit frequency per unit sample area per unit solid angle) is given by Planck's formula for blackbody emission S_{bb}(ν, T) = (2*h*ν^{3}/*c*^{2})/[exp(*h*ν/k_{B}T) − 1] multiplied by the nanoparticle absorption efficiency. Importantly, while this result is derived under the assumption of thermodynamic equilibrium, it can also be applied under a broad range of nonequilibrium conditions.^{36} The absorption efficiency η(ν) is defined as the ratio between absorption cross section and physical surface area, and can be well above one at frequencies near a nanoparticle absorption resonance. For the devices under study, it can be written as η_{pl}(ν) = α_{pl}(ν)S_{tot}/S_{gr}, where S_{gr} is the graphene surface area, S_{tot} is the total device area (including both graphene nanoribbons and the empty spaces in between), and α_{pl}(ν) is the fraction of incident optical power that is absorbed by the plasmonic oscillations.

Based on this model, the plasmonic emission spectrum and output power can be computed as follows. First, FDTD simulations are carried out to calculate the graphene plasmonic absorption α_{pl}(ν) for an input plane wave of variable polar and azimuthal angles of incidence (θ and ϕ, respectively) and different states of polarization (e.g., *s* and *p*). Next, the output power density per unit frequency per unit sample area ρ_{out}(ν, T) is obtained by integrating the resulting spectral radiance over all angles, i.e.,

where the brackets ⟨⟩ indicate averaging over two orthogonal incident polarizations. With this formulation, the wall-plug efficiency ΔP_{out}/P_{in} can finally be expressed as

where *p*_{in} = P_{in}/S_{gr} is the input electrical power density, T_{base} is the sample temperature for P_{in} = 0, and T_{el} is the elevated temperature of the graphene electron gas under bias. With *p*_{in} given, the latter parameter can be evaluated using the heat-transfer model of Refs. 37–39, where “supercollisions” involving disorder-assisted acoustic-phonon scattering are identified as the dominant cooling mechanism for hot carriers in graphene. In this model, $pin=A[Tel3\u2212Tlatt3]$, where A is the supercollision coupling constant (proportional to carrier density N^{37}). Furthermore, the lattice temperature under bias is computed as T_{latt} = T_{base} + (r_{G–OX} + t_{OX}/κ_{OX})*p*_{in}, where r_{G–OX} = 4 × 10^{−8} K m^{2}/W is the interface thermal resistance between graphene and the SiO_{2} gate oxide,^{40} t_{OX} = 300 nm is the oxide thickness, and κ_{OX} = 1.3 W/m/K is the thermal conductivity of SiO_{2}. This expression for T_{latt} neglects the lateral spreading of the heat flow in the supporting materials and the contribution from the Si substrate. Nevertheless, its predictions are in very good agreement with more detailed calculations and measurements of the lattice temperature in the same graphene FET device structure.^{40}

The results of this analysis for the experimental device of Fig. 1(b) are summarized in Fig. 2. Figure 2(a) shows the calculated absorption spectrum α(ν) of this device illuminated at normal incidence with light polarized parallel (*y*) and perpendicular (*x*) to the nanoribbons. In these FDTD simulations, the sample quality is quantified by the carrier scattering lifetime τ used in the graphene conductivity model or, equivalently, by the carrier mobility $\mu =qvF2\tau /EF$, where v_{F} ≈ 1 × 10^{8} cm/s and E_{F} are the graphene Fermi velocity and Fermi energy, respectively. The spectra plotted in Fig. 2(a) were computed for μ = 2000 cm^{2}/V s, which is appropriate for typical CVD graphene samples deposited on SiO_{2} using PMMA.^{41} The shapes of these spectra reflect the well-established polarization properties of GPPs in nanoribbons,^{7,10,24,34,35} where only the perpendicular component of the incident electric field can experience resonance absorption by coupling to the fundamental plasmonic oscillations. In contrast, a weaker absorption coefficient decreasing monotonically with increasing frequency is observed for the parallel field component, indicative of free-carrier absorption. As a result, in the nanoribbon geometry, the polarization averaging expressed by the brackets in Eq. (1) is equivalent to multiplication of the *x*-polarized absorption by a factor of 1/2. It follows from Fig. 2(a) that the polarization-averaged peak absorption ⟨α_{pl}(ν_{pl})⟩ of the device under study at normal incidence is ∼0.025, corresponding to a rather small absorption efficiency ⟨η_{pl}(ν_{pl})⟩ = ⟨α_{pl}(ν_{pl})⟩S_{tot}/S_{gr} = 0.047.

In Fig. 2(b), the polarization-averaged peak absorption efficiency ⟨η_{pl}(ν_{pl})⟩ is plotted as a function of the polar angle of incidence θ on the *x–z* and *y–z* planes [relative to the system of coordinates defined in Fig. 2(a)]. A monotonic decrease in absorption with increasing illumination angle is observed on both planes, mostly arising from the corresponding increase in the projected area of the incident wave on the sample plane as well as polarization effects. The full angular dependence of ⟨η_{pl}(ν_{pl})⟩ on both θ and ϕ is displayed in Fig. 2(c). For computational convenience, this color map was calculated with a reciprocity-based method, whereby an electric dipole source at ν_{pl} is positioned on a nanoribbon and its far-field radiation pattern is computed in the air above. By reciprocity,^{42} this pattern is proportional to the local field intensity (and therefore the absorption) produced at the dipole location by an incident plane wave as a function of illumination angles. The proportionality factor between the calculated radiation pattern and the absorption efficiency ⟨η_{pl}(ν_{pl})⟩ is then determined by comparing the horizontal and vertical line cuts of Fig. 2(c) with the traces of Fig. 2(b). Similar results for the angular dependence of ⟨η_{pl}(ν)⟩ were obtained at all other frequencies ν across the plasmonic absorption spectrum.

Using these simulation results in Eq. (1), we can finally calculate the plasmonic emission spectrum Δρ_{out}(ν) = ρ_{out}(ν, T_{el}) − ρ_{out}(ν, T_{base}) of the device of Fig. 1(b) under the same experimental conditions (i.e., T_{base} = 80 K, P_{in} = 0.2 W, and N = 5.4 × 10^{12} cm^{−2}). The supercollision coupling constant A used to evaluate the electronic temperature under bias T_{el} is computed with the same material parameters described in Ref. 24. A value of T_{el} = 196 K is obtained in this case, whereas the corresponding lattice temperature T_{latt} is 83 K. The resulting emission spectrum, plotted in Fig. 2(d), is in reasonably good agreement with the experimental data of Fig. 1(b). Specifically, the calculation results show a peak emission frequency ν_{pl} = 8.0 THz, a spectral linewidth δν = 2.8 THz (full width at half maximum), and a peak power density Δρ_{out}(ν_{pl}) = 14.5 *μ*W/THz/cm^{2}, whereas the measured values are ν_{pl} = 7.9 THz, δν = 2.0 THz, and Δρ_{out}(ν_{pl}) = 9.1 *μ*W/THz/cm^{2}. Altogether, this comparison is quite satisfactory, particularly given the absence of freely adjustable fitting parameters in the simulations, and therefore it substantiates the validity of the underlying theoretical model.

In order to elucidate the key limiting factors of the device structure considered so far, in Fig. 3, we plot its normal-incidence polarization-averaged peak absorption ⟨α_{pl}(ν_{pl})⟩ and corresponding peak absorption efficiency ⟨η_{pl}(ν_{pl})⟩ as a function of carrier mobility μ (bottom axis) and carrier scattering rate 1/τ (top axis). The shape of this trace can be explained with a general argument based on the coupled-mode theory:^{43,44} the coupling between external radiation and the plasmonic resonance is maximum when the radiative decay rate of the plasmonic oscillations Γ_{rad} is equal to their nonradiative damping rate Γ_{nr}. Away from this critical coupling condition, the absorption decreases regardless of whether Γ_{rad} < Γ_{nr} (under-coupling) or Γ_{rad} > Γ_{nr} (over-coupling). In a 2D electron gas, nonradiative damping of the plasmonic oscillations involves the same scattering mechanisms that also limit the mobility (impurities, defects, phonons, and other carriers). Therefore, Γ_{nr} can be taken to be equal to 1/τ, which is inversely proportional to μ as discussed above. The radiative decay rate Γ_{rad} is instead mostly limited by the large wavelength mismatch between the graphene surface plasmons and free-space radiation. The simulation results of Fig. 3 indicate that in the present device structure, Γ_{rad} is ∼3.4 × 10^{12} s^{−1}, corresponding to an optimal mobility for critical coupling μ* ≈ 11 000 cm^{2}/V s. In the experimental sample of Fig. 1(b), the actual mobility is significantly smaller (about 2000 cm^{2}/V s near 80 K), and therefore its plasmonic emission is strongly under-coupled.

These results suggest three key design guidelines to increase the emission efficiency and correspondingly enable operation at higher temperatures. First, optimal device structures should be designed to operate as close as possible to critical coupling. In principle, this goal could be achieved by employing ultrahigh-quality graphene sheets embedded in hexagonal-BN^{45,46} or partially suspended,^{47,48} where room-temperature mobilities approaching the optimal value μ* of Fig. 3 are possible. However, at present, these samples tend to be quite small and cannot be readily scaled to the mm^{2} dimensions required to obtain significant output radiation power. Alternatively, the graphene nanostructures may be combined with additional photonic and/or plasmonic antennas that can effectively “funnel” incident radiation in and out of the GPPs,^{19,49} in order to increase their radiative decay rate Γ_{rad} and thus lower the optimal mobility μ*. Second, the graphene dielectric environment can be engineered to increase the peak plasmonic absorption at critical coupling (α* ≈ 0.046 in Fig. 3) toward its maximum possible value of 100%. Aside from the polarization properties described above, α* here is limited by the presence of both transmission and reflection channels for the incident light, in addition to absorption, so that perfect absorption may be achieved by suppressing these alternative processes.^{50,51} Finally, according to Eq. (1), the output radiation power can be further increased by maximizing the geometrical ratio S_{tot}/S_{gr} for fixed absorption ⟨α_{pl}(ν_{pl})⟩, i.e., by scaling the sample area occupied by each nanoribbon to match its absorption cross section.

## III. CRITICALLY COUPLED GRAPHENE PLASMONIC OSCILLATORS

Figure 4 shows a general device geometry designed following the guidelines just presented, where each graphene nanoribbon is surrounded symmetrically by two THz antennas consisting of parallel Au rectangular patches. These metallic antennas support resonant modes at free-space wavelengths λ proportional to their lateral dimension L, with a moderate degree of spatial confinement λ/L ∼ 10^{52} (as opposed to λ/*w* ∼ 100 for the graphene nanoribbons). Therefore, they can effectively mediate the wavelength mismatch limiting the GPP radiative decay rate,^{49} without at the same time adding any significant absorption losses by virtue of the small penetration depth of THz light in metals. For fixed lateral dimensions *w* and L, the coupling between each nanoribbon and its neighboring antennas is controlled by their separation *s*, whereas the geometrical ratio S_{tot}/S_{gr} can be optimized by selecting the array period P. Additionally, an open vertical cavity configuration is also employed in the device of Fig. 4, consisting of a Au back mirror, a low-doped Si layer of thickness *h* (which can also serve as the back conductor of the gate capacitor), and a thin SiO_{2} gate dielectric. Under external illumination, the back mirror can block the transmission of the incident light and at the same time suppress reflection through interference (with the proper choice of graphene-mirror separation), and therefore maximize the absorption.^{50,51} A similar vertical cavity configuration was originally developed in the context of radio-wave engineering, where it is referred to as the Salisbury screen.^{53}

Design optimization for the device structure of Fig. 4 involves several geometrical parameters, which collectively determine the absorption efficiency. In the work presented below, we have fixed the Au-antenna and SiO_{2}-gate-dielectric thicknesses to 100 nm and 300 nm, respectively, and varied the remaining parameters in order to maximize the plasmonic emission near 8 THz (for a direct comparison with the reference device of Figs. 1–3). Representative design simulation results are plotted in Figs. 5 and 6. Figure 5(a) shows the resonance frequency of minimum reflection for the coupled Au-antenna/vertical-cavity system (without the graphene nanoribbons) under *x*-polarized normal-incidence illumination, for different values of the Si thickness *h* and antenna width L. The antennas in these simulations are arranged periodically on the SiO_{2} surface with a fixed edge-to-edge separation of 600 nm. In this configuration, reflection is partially suppressed via destructive interference between the waves reflected at the top and bottom surfaces of the vertical cavity (i.e., at any frequency where the total roundtrip phase shift through the cavity is π), combined with absorption inside the device (mostly in the SiO_{2} layer as discussed below). As shown in Fig. 5(a), the resulting resonance frequency decreases with increasing *h* and L, due to the resulting variations in the propagation phase through the Si layer and in the scattering phase of the Au antennas, respectively. Based on these results, we select *h* = 1.5 *μ*m and L = 3.32 *μ*m, which produces a pronounced reflection dip at the target operation frequency of 8 THz [inset of Fig. 5(a)].

Next, we introduce the graphene nanoribbons near the Au patch antennas and evaluate their plasmonic absorption spectra for different values of the graphene Fermi energy E_{F} and mobility μ. Initially, we consider a variation of the device geometry of Fig. 4 where a single antenna is located symmetrically between consecutive nanoribbons [inset of Fig. 5(b)], with 400-nm ribbon width *w* and 100-nm ribbon/antenna separation *s*. Figure 5(b) shows the resulting plasmonic absorption α_{pl} under *x*-polarized normal-incidence illumination, plotted as a function of frequency ν and Fermi energy E_{F} for μ = 2000 cm^{2}/V s. The shape of this color map reflects the coupling between the nanoribbon fundamental plasmonic excitation, whose frequency increases with the square root of E_{F}^{7,10,24,34,35} (red dashed curve), and the antenna resonance which is independent of E_{F} (orange dashed line). To further elucidate the absorption properties of this device, in Fig. 5(c), we plot the different contributions to its absorption spectrum for constant E_{F} = 0.23 eV. As shown by these simulation results, the main competing absorption mechanism limiting the peak value of α_{pl}(ν) is provided by the SiO_{2} gate dielectric, which incidentally is also responsible for an absorption peak near 14 THz associated with a SiO_{2} surface-phonon resonance.^{10} At the same time, free-carrier absorption in the Si film is found to be quite small in the spectral region of interest, even though the Si doping level used in these simulations (1 × 10^{16} cm^{−3}) is already large enough for back-gating purposes.^{24}

In any case, the key conclusion that emerges from Fig. 5(c) is that the device under study can provide near-complete absorption (for *x*-polarized, normally incident light) at the GPP resonance frequency ν_{pl}, despite its relatively low graphene mobility of 2000 cm^{2}/V s. In fact, this value essentially coincides with the optimal mobility μ* for critical coupling of the same device, as shown in Fig. 5(d) where the polarization-averaged peak absorption ⟨α_{pl}(ν_{pl})⟩ and corresponding peak absorption efficiency ⟨η_{pl}(ν_{pl})⟩ are plotted as a function of μ and 1/τ. As a result of the combined action of the vertical cavity and THz antennas, we find that the GPP radiative decay rate Γ_{rad} is increased from 3.4 × 10^{12} s^{−1} in the “bare” nanoribbons of Fig. 3 to nearly 20 × 10^{12} s^{−1} in the present device, leading to a proportional decrease in μ* to about 2000 cm^{2}/V s. Importantly, this value of mobility is readily accessible even at room temperature with large-area CVD-grown graphene transferred on Si/SiO_{2} substrates.^{41} Therefore, maximum plasmonic absorption can be achieved with this device configuration without the need for ultrahigh-mobility graphene, only limited by the nanoribbon polarization properties and the above-mentioned competing processes in the SiO_{2} gate dielectric. Furthermore, the absorption efficiency ⟨η_{pl}(ν_{pl})⟩ at μ = 2000 cm^{2}/V s is as large as 3.8 (vs 0.047 for the reference device of Figs. 1–3) and remains well above unity for mobilities as small as a few 100 cm^{2}/V s. As a result, significantly stronger plasmon-enhanced thermal emission can be expected from the present device, even under ambient conditions at room temperature.

The absorption efficiency can be further enhanced by increasing the center-to-center inter-nanoribbon separation P in order to optimize the geometrical ratio S_{tot}/S_{gr}. To this end, we introduce two separate antennas between neighboring nanoribbons (as in Fig. 4) with fixed separation *s* = 100 nm between each ribbon and its surrounding antennas, so as to preserve their strong coupling responsible for the large decrease in μ* just discussed. Figure 6(a) shows the normal-incidence absorption ∫dν⟨α_{pl}(ν)⟩ and absorption efficiency ∫dν⟨η_{pl}(ν)⟩ integrated over all frequencies and plotted as a function of the inter-nanoribbon separation P. All other geometrical parameters used in these simulations are the same as in Fig. 5(d), with μ = 2000 cm^{2}/V s and E_{F} = 0.27 eV (corresponding to a carrier density N = 5.4 × 10^{12} cm^{−2}, as in the reference device of Figs. 1–3). As expected, the absorption decreases monotonically with increasing P, due to the corresponding decrease in the fraction of incident power intercepted by the nanoribbons. At the same time, however, the absorption efficiency initially increases with P, indicating that the absorption cross section of each nanoribbon along the *x* direction is initially larger than the separation between consecutive ribbons. Based on these simulation results, we can maximize the absorption (and therefore the radiation) efficiency by setting P = 20 *μ*m, which is approximately half the free-space wavelength at the peak emission frequency of 8 THz.

In order to evaluate the plasmonic emission spectrum and output power of the resulting device using Eqs. (1) and (2), we finally compute ⟨η_{pl}(ν)⟩ as a function of direction of propagation of the incident light. Representative results for different frequencies ν are shown in Figs. 6(b)–6(d), where once again the color maps of absorption efficiency vs polar and azimuthal illumination angles [panels (c) and (d)] were obtained with the reciprocity-based method described above. Compared to the reference device of Fig. 2, more complex patterns are observed in these plots, with narrower peaks and in some instances sidelobes depending on the frequency. This behavior originates from interference at the location of each nanoribbon between the incident (or radiated) light and the light that is correspondingly scattered by the neighboring ribbons. More specifically, in the 3D simulations of Figs. 6(c) and 6(d), we consider a small portion of the entire device with only five adjacent repeat units, so that the largest separation between different nanoribbons is 4 × 20 = 80 *μ*m. This value is smaller than the estimated coherence length of the device plasmonic emission (127 *μ*m in free space, as discussed below), so that the interference effects observed in these plots are consistent with the partial spatial coherence of the emitted light. In passing, it should be noted that the use of periodically patterned surfaces to produce directional and collimated thermal emission is already well established.^{54}

The plasmonic emission spectrum Δρ_{out}(ν) = ρ_{out}(ν, T_{el}) − ρ_{out}(ν, T_{base}) computed with Eq. (1) using the absorption data just described is shown in Fig. 7(a) for P_{in} = 100 W and S_{gr} = 1 mm^{2} (importantly, the ability of graphene to withstand such relatively large electrical power densities has already been reported^{29,31,40}). The base temperature T_{base} used in these calculations is 300 K, with the corresponding lattice and electronic temperatures under bias estimated at T_{latt} = 327 K and T_{el} = 448 K, respectively. A relatively sharp emission spectrum is observed in Fig. 7(a) peaked at 8.5 THz, near the “bare” plasmonic resonance of the nanoribbons, with a somewhat asymmetric line shape resulting from the coupling with the modes of the adjacent Au antennas. The emission linewidth δν is 2.1 THz, corresponding to a coherence length *l*_{c} = (4ln2/π)(*c*/δν)^{55} of 127 *μ*m as mentioned above. The total output power ΔP_{out} in Fig. 7(a) is 0.28 mW for a wall-plug efficiency ΔP_{out}/P_{in} of 2.8 × 10^{−6}. For comparison, using the same procedure with the reference bare nanoribbon device of Fig. 2 under the same conditions (T_{base} = 300 K, *μ* = 2000 cm^{2}/V s, N = 5.4 × 10^{12} cm^{−2}, P_{in} = 100 W, and S_{gr} = 1 mm^{2}), we find a much lower output power ΔP_{out} of 1.4 *μ*W. We can therefore conclude that the use of structures designed to enable critical coupling and perfect absorption can increase the plasmonic radiation efficiency by two orders of magnitude. Finally, our simulation results also show that for fixed P_{in}, the output power increases with increasing graphene area S_{gr} [Fig. 7(b)], despite the corresponding decrease in electronic temperature T_{el} [Fig. 7(c)]. In this respect, it should be noted that CVD graphene can be synthesized in extremely large, centimeter-scale films,^{56} although the active area of the present devices is also limited by their large geometrical ratio S_{tot}/S_{gr}.

## IV. CONCLUSION

We have presented a detailed theoretical model of plasmon-enhanced graphene thermal emission based on the recent demonstration of current-driven THz radiation from graphene nanoribbons.^{24} The same model has then been applied to design a more complex device structure with significantly enhanced radiative efficiency. This design leverages the distinct ability of resonant nanostructures to provide absorption cross sections much larger than their physical area (especially under conditions of critical coupling), combined with the ultrastrong light–matter interactions enabled by graphene plasmonic excitations. The same general framework could also be further expanded for additional performance improvements. Specifically, the devices considered in the present study are based on relatively simple plasmonic oscillators (nanoribbons) consisting of standard CVD-grown graphene on Si/SiO_{2} substrates, with readily accessible carrier densities of a few 10^{12} cm^{−2}. A larger absorption (up to twofold) could be achieved with graphene nanostructures featuring a polarization-independent plasmonic response, for example, circular or slightly elliptical disks connected to one another by narrow ribbons to enable current flow and gating.^{57} Plasmonic absorption in graphene can also be enhanced by increasing the carrier density, e.g., using a thin ion-gel top-gate architecture as in the initial measurement of THz GPPs in nanoribbons,^{7} where densities exceeding 15 × 10^{12} cm^{−2} were obtained. The use of double-layer graphene is similarly advantageous to increase the accessible doping range.^{58,59} In addition, samples combining large areas with higher room-temperature mobilities (e.g., ∼10 000 cm^{2}/V s for CVD-grown graphene transferred on h-BN^{60}) could also be envisioned as a way to simultaneously increase the peak power density and decrease the spectral linewidth of the output light. This approach is particularly interesting as a way to enable sensitive gate tuning of the radiation spectra, which is currently limited by the resonant nature of the Au-antenna/vertical-cavity system combined with the broader plasmonic emission in lower-mobility samples [e.g., see Fig. 5(b)]. A narrower spectral linewidth would also increase the spatial coherence of the output light and, therefore, could be exploited to enhance directionality of the radiation pattern.

Altogether, these results indicate that technologically significant THz output power levels approaching the mW range can be achieved with properly designed graphene plasmonic oscillators under continuous-wave room-temperature operation. The required input electrical powers (on the order of 100 W) are large but accessible and compatible with the high-current-carrying capabilities of graphene.^{29,31,40} It should also be noted that current solid-state technologies for THz light emission, most notably quantum cascade lasers (QCLs),^{61} generally require cryogenic cooling, which significantly limits their portability and integration in practical systems. More recent solutions, including thermoelectrically cooled THz QCLs^{62,63} and intracavity difference-frequency generation in room-temperature mid-infrared QCLs,^{64} are instead mostly limited to pulsed emission with wall-plug efficiencies smaller than or comparable to the graphene devices described above. These devices, therefore, appear to be promising for a wide range of practical applications in THz imaging and spectroscopic sensing. Specific advantages include direct compatibility with Si microelectronics, relatively simple manufacturing (not requiring any complex heteroepitaxial growth), and the ability to operate across the entire THz spectrum, including frequencies that are fundamentally inaccessible to III–V semiconductor materials due to reststrahlen absorption (e.g., ∼7–9 THz in As-based materials).

## ACKNOWLEDGMENTS

The FDTD simulations were performed using the Shared Computing Cluster facility at Boston University. Y.L. acknowledges partial support from the National Science Foundation under Grant No. ECCS-1711156.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: SIMULATION METHODS

All the FDTD simulations reported in this work were performed using the Lumerical FDTD Solutions software package. In these simulations, graphene is treated as a 2D interface with a complex frequency-dependent conductivity based on the semi-classical Drude model. All other relevant materials (Si, SiO_{2}, and Au) are described with their complex dielectric functions. Specifically, for Si and SiO_{2}, we used the built-in database of the FDTD software; for Au, we used data from Ref. 65, which cover a broad spectral range including the THz frequencies of interest in this work. The absorption spectra are computed via 2D simulations, where the device under study is illuminated with a linearly polarized plane wave incident from the air above. Periodic boundary conditions are applied to the lateral boundaries of the simulation region, whereas perfectly matched layers (PMLs) are used in the top and bottom surfaces. The individual contribution to the total absorption from the graphene nanoribbons, or any other layer of interest, is isolated using a built-in function of the FDTD solver. A broadband input light pulse covering the spectral range of interest is used in all simulations, except for those involving incidence at oblique angles [as in Figs. 2(b) and 6(b)] where the use of a single-frequency source is required. Finally, the far-field radiation patterns are computed via 3D simulations, where the computational domain consists of five adjacent repeat units along the *x* direction with sufficiently large extent (20 *μ*m) along the *y* direction and PMLs on all boundaries. In these simulations, the emitted light is provided by an electric dipole oriented along the *x* direction, located on the graphene nanoribbon at the center of the simulation region. The radiation patterns of *y*- and *z*-oriented dipole sources were also computed but found to be negligibly small by comparison, consistent with the polarization properties of the plasmonic resonances under study.