The device we used to amplify THz waves is shown in Fig. 1A. Our device is composed of a graphene layer interposed between 2 dielectric media. The top dielectric layer is overlaid by a metal gate. The graphene layer supports surface plasmons. The attractiveness of employing graphene to support surface plasmons lies in the fact that the carrier concentration within graphene can be continuously tuned from nearly zero to as high as 10¹³ cm⁻² through the application of an external voltage. This important ambipolar doping effect has unveiled the enticing possibility of achieving real-time control of surface plasmons [54].

In the proposed mechanism, the presence of a drift velocity induces a strong asymmetry in the electron–phonon scattering rate across the Fermi surface. This asymmetry leads to electron redistribution in phase space during transport, enhancing local charge concentration perturbations and thereby effectively amplifying the plasmons. To better illustrate our findings, we physically model the system as described below.

Two types of electromagnetic surface waves, namely, the transverse electric (TE) and transverse magnetic (TM) surface modes, are capable of propagating in graphene. Arising from longitudinal charge concentration oscillations, graphene plasmons exhibit a longitudinal electric field component, classifying them as TM modes. In the THz and far-infrared regimes, graphene’s conductivity is dominated by intraband carrier dynamics, which inherently facilitates TM-mode propagation by enabling efficient coupling between electromagnetic fields and charge concentration fluctuations [55]. Here, the Boltzmann transport equation [56–58] is adopted to describe the motion of electrons in phase space:

Here, v_F is the Fermi velocity, k = (k_x, k_y) is the wavevector, e is the electron charge, E_x denotes the longitudinal electric field, ħ is the reduced Planck’s constant, f(x,k,t) is the distribution function of electrons, and S_e-ph is the collision integral that accounts for electron–phonon scattering. Combined with Maxwell’s equations, the graphene plasmon dispersion is derived:

Here, g = g_sg_v stands for the degeneracy of electrons, g_s = 2 is spin degeneracy and g_v = 2 is the valley degeneracy, μ is chemical potential, and ε₀ denotes the permittivity of free space. {\gamma }_1=\sqrt{\left(q^2-{\omega }^2{\epsilon }_1/c^2\right)}is the constant of evanescent decay of the surface wave field in the transverse direction from graphene; d₁ is the height of the top media; q and ω are the wavevector and angular frequency of surface plasmon, respectively; k_F is the Fermi wavevector; Г is the scattering rate between electron and phonon; and β = v_x/v_F is the relative velocity. The details of theoretical derivation are provided in Methods. By solving the dispersion relation [Eq. (2)] numerically, the complex oscillation frequency of graphene plasmon \omega ={\omega }^{\prime }+{\mathit{i\omega }}^{″} is obtained. Its real and imaginary parts represent the resonances and increments (amplification efficiency).

Figure 1A illustrates the THz amplification due to graphene plasmon instability induced by drift velocity in a single-layer gated graphene system. In the proposed mechanism, plasmons propagating along the drift velocity (downstream plasmons) attenuate gradually, while plasmons propagating against the drift velocity (upstream plasmons) amplify. Figure 1B shows that the normalized electric field of plasmons is sinusoidal along the x direction, with its maximum amplitude occurring in the graphene layer (z = 0) and decaying in the medium. This reflects the strong confinement of the graphene plasmons. Figure 1C depicts the energy diagrams of graphene without and with electron drift velocity. It can be observed that the Fermi energy is constant without drift velocity. However, the Fermi surface tilts within the Dirac cone when the electrons have a drift velocity, causing the Fermi energy to change with momentum. From Fig. 1D, we observe that the introduction of drift velocity markedly alters the plasmon dispersion, leading to nonreciprocal plasmon propagation that has been experimentally verified [30,31]. Due to the Doppler effect, the downstream plasmon exhibits a blueshift at a fixed wavevector accompanied by an enhanced group velocity. Conversely, the upstream plasmon is redshifted with a reduced group velocity. Furthermore, the results reveal that the frequency shifts become more notable as the drift velocity increases, indicating a direct correlation between the Doppler effect and the frequency of the amplified mode. As shown in Fig. 1E, positive increments occur for upstream plasmon propagation, indicating the amplification in the THz range. The amplification efficiency can be modulated by the electron drift velocity. As the drift velocity gradually increases from 0.2v_F to 0.8v_F, the amplification efficiency increases from a negative value to 0.16 ps⁻¹. Further, the positive increment persists over a broad range of wavevectors, suggesting the characteristic of ultra-broadband amplification. The maximum increment occurs in the long wavelength region, corresponding to several THz. This may imply that the plasmons with several THz will dominate in the amplification process. It is expected to provide a possibility for the scarce THz sources (0.1 to 10 THz). Figure 1F shows that electron concentration is another key factor that can markedly regulate the amplification efficiency. As the electron concentration increases from 1 × 10¹² cm⁻² to 5 × 10¹² cm⁻², the amplification efficiency increases from about 0.15 ps⁻¹ to close to 1 ps⁻¹, which means a decrease in propagation distance for THz amplification. We also observe the peak shifts to larger wavevectors. This shift may result in the plasmons with higher frequency becoming dominant.

To reveal the properties of the instability, we investigate the threshold required for amplification. This analysis explores a broad range of electron concentrations, encompassing typical doping levels for graphene (0.1 to 0.5 eV), and considers drift velocities ranging from zero to the Fermi velocity v_F. Figure 2A shows the region where THz amplification is achieved under the combined influence of electron concentration and electron drift velocity. It indicates that the drift velocity threshold for plasmon instability initially decreases and then increases with increasing electron concentration. The minimum threshold occurs at approximately n = 2 × 10¹² cm⁻², indicating that only a relatively low drift velocity is needed for amplification at this concentration. This is particularly advantageous for practical applications, as it relaxes the limitation for high electron mobility. To explain this phenomenon, we choose 3 dashed lines in Fig. 2A and plot their maximum scattering rate differences (ΔГ = Г_max − Г_min) and the relative differences ΔГ/Г_max on the Fermi surface. The results are shown in Fig. 2B. In the THz amplification region, the scattering rate differences ΔГ are all on the order of 1 and 10 ps⁻¹. Meanwhile, the relative differences ΔГ/Г_max are also large, approximately greater than 0.8. This indicates that the key to achieving THz amplification is to create large relative scattering rate differences on the Fermi surface. Further, we investigate the scattering rates of electrons with acoustic and optical phonons. Figure 2C shows that the strength of the electron–optical phonon interactions is much higher than that of the electron–acoustic phonon interactions. Combined with Fig. 2B, we can conclude that the THz amplification is dominated by the electron–optical phonon interactions. We also demonstrate that although the instability is determined jointly by 3 optical phonon modes, the contribution from the phonons near the K point dominates due to its higher scattering rate with electrons (see detailed analysis in the Supplementary Materials).

Figure 2D tracks the amplitude evolution of the electron density perturbation over one period, showing an initial decay followed by a sustained growth. To elucidate the amplification mechanism, we examine the evolution of electron momentum-space distribution at the wave crest, as presented in Fig. 2E and F. The results indicate a growing electron clustering in low-scattering-rate regions near the Fermi surface at negative momenta (−k_x), leading to the enhancement of the wave amplitude. Based on these findings, we propose the following explanation. Initially, the drift velocity induces a nonuniform scattering rate distribution across the tilted Fermi surface. Perturbations established in high-scattering-rate regions are rapidly damped, causing the initial wave attenuation. Meanwhile, the electrons are redistributed in the momentum space and accumulate preferentially in low-scattering-rate regions. Although the individual electron velocity is lower than the phase velocity of the wave, the collective motion of the electron cluster synchronizes with the wave crest, thereby continuously delivering energy to the wave. Once the rate of electron cluster formation exceeds the wave-damping rate (i.e., the threshold in Fig. 2A is surpassed), net amplification occurs. The excited plasmon must therefore acquire negative momentum from the clustered electrons through momentum conservation, confirming its upstream nature.

Based on the above physical explanation, we believe that graphene has 2 unique advantages over other plasmonic materials or 2D systems. Specifically, graphene’s linear energy dispersion enables the formation of electron flows moving against the drift, especially under conditions of high drift velocities. This counterflow plays a crucial role in making the amplification of upstream plasmons possible. Additionally, the linear energy dispersion leads to a broad carrier distribution in momentum space, resulting in a large variation in energy at the tilted Fermi surface. This characteristic amplifies the variations in scattering rates across the Fermi surface, thereby enhancing the proposed amplification mechanism.

Compared to the recent dissipative instability in single-layer graphene [43], the instability proposed in this work similarly employs electron drift velocity and momentum relaxation to achieve non-Cherenkov amplification, with similar tunable increments (0.1 to 1 ps⁻¹) in the THz range. However, their results differ due to the distinct physical models. The dissipative instability amplifies downstream plasmons with phase velocities near the electron drift velocity, whereas the present mechanism amplifies upstream plasmons and exhibits pronounced sensitivity to momentum-space variations in the relaxation rate. The different theoretical approaches also determine their respective operational conditions. The dissipative instability is valid in the hydrodynamic regime. It requires that the mean free path of electron–electron scattering l_ee becomes much shorter than that of the electron–phonon scattering and the system dimensions [59]. In contrast, the electron–phonon scattering dominates over the electron–electron scattering in the present mechanism, which means that the mean free path of electron–phonon scattering l_e-ph is shorter than l_ee. This difference implies that l_ee markedly affects the applicability of the 2 mechanisms. According to Ref. [60], l_ee in single-layer graphene varies with carrier density n and temperature T as l_ee ∝ n^1/2/T². Therefore, in practical high-mobility graphene field-effect transistors, it is possible to switch between the 2 mechanisms by adjusting the gate voltage and temperature. Furthermore, as will be discussed in the subsequent section, the nonlinear dynamics associated with the proposed mechanism offer a promising route toward THz frequency comb generation, potentially extending its utility to high-precision spectroscopy and broadband communications.

To study the amplification process of graphene plasmon, in particular the influence of nonlinear effects on THz amplification, we take advantage of the semi-Lagrangian method [53] to solve the Boltzmann transport equation [Eq. (1)]. The semi-Lagrangian method solves kinetic equations by backwardly tracing particle trajectories in phase space on a fixed grid. It employs interpolation with time-splitting technique to update the distribution function at each time step, avoiding grid distortion while retaining stability for calculation. The detailed procedures of the semi-Lagrangian method are given in the Supplementary Materials.

Here, we briefly describe the construction of the grids. We set Δx ≈ 5 nm, Δk_x ≈ 3 μm⁻¹, and Δk_y ≈ 5 μm⁻¹, where Δx, Δk_x, and Δk_y are the spatial grid steps along x, k_x, and k_y directions, respectively. Given that the electron drift velocity is aligned with the x direction, the electron distribution function exhibits enhanced broadening in momentum space along k_x. To resolve this anisotropic feature, we employ a smaller grid spacing (Δk_x) in the k_x direction. A nonuniform temporal step size, Δt ϵ [0.1,5] fs, is implemented, as the nonlinear stage evolution introduces higher harmonics, which require shorter time steps for accurate simulation. Additionally, periodic boundary conditions are imposed on the boundaries along the x direction. Multiple simulation results show that the above parameter settings maintain the convergence and stability of the solution. In the amplification process, the longitudinal electric field of graphene plasmons has the following form:E_x=E_{x-\mathit{\text{initial}}}{e}^{{\omega }^{\prime \prime }t}{e}^{-i{\omega }^{\prime }t}. Here, E_x-initial is the amplitude of the longitudinal electric field at the initial time. It is determined by the applied perturbation and can be calculated by Eq. (7). It is worth mentioning that the saturated electric field strength is independent of the initial perturbation, which is mainly determined by the equilibrium electron concentration and electron drift velocity. (The detailed analyses are provided in Section S4.) So, it is appropriate to introduce a dimensionless parameter ζ = ln (E_x/E_x-initial) as the relative amplitude of the electric field to characterize the linear amplification process.

In Fig. 3A and B, we introduce a sinusoidal perturbation to the electron concentration whose amplitude is 0.01% of the equilibrium concentration. The results show that the amplification process of graphene plasmons is consistent at different electron concentrations. The relative amplitude of the longitudinal electric field first begins to amplify with a specific increment in accordance with the solution of dispersion and then becomes saturated. The differences brought about by different electron concentrations are mainly reflected in the following 2 aspects. The first is the time required for amplification to saturation, which is about 12.5 and 2.5 ps, respectively. The second is the peak and gain of the electric field of the plasmons. The electric field can peak at 10⁶ and 10⁷ V/m, and the gain can reach about 58 and 70 dB, respectively. Therefore, we believe that high electron concentrations are more favorable for ultrafast and high-gain THz amplification. To capture the changes in the spectrum more clearly, we choose the case where the electron concentration is relatively low, i.e., n = 1 × 10¹² cm⁻², because in this condition the amplification process is slower and the spectrum does not change too rapidly. Figure 3C shows that in the linear stage, the plasmon frequency is primarily determined by the wavevector of the initial perturbation. During this stage, the plasmon can be approximated as a single-frequency wave, with its amplitude gradually increasing at a constant rate. This mode corresponds to the fundamental mode. When the fundamental mode is amplified adequately, as shown in Fig. 3D and E, higher harmonics are excited, signaling the transition from the linear stage to the nonlinear stage. The nonlinear effect originates from the third term on the left-hand side of Eq. (1), which contains the products of any-order perturbations of the electric field and electron distribution function. These higher harmonics align with the frequency-doubling line, indicating that their frequencies and wavevectors are integer multiples of the fundamental mode. At this time, amplification of the higher harmonics proceeds simultaneously. Notably, these higher-order modes generate a THz frequency comb with a broad spectrum, extending up to 25 THz. Subsequently, these higher harmonics, which initially follow the frequency-doubling line, gradually approach the dispersion, as shown in Fig. 3F. During this process, these frequencies decrease, causing the frequency comb to shift toward the lower frequency band. Meanwhile, these plasmon modes compete with each other and eventually evolve to the dominance of the higher harmonics. In general nonlinearities, the intensity of the higher harmonics should be much less than that of the fundamental mode due to the low conversion rate. However, during the present amplification process, the higher harmonics are gradually enhanced and eventually dominate. In this regard, we infer that this result is due to the property of the ultra-broadband spectrum of the plasmon instability. This means that all the higher harmonics falling within the frequency region of the instability can be amplified.

Additionally, we have investigated the influence of static electric fields E₀ along the −x direction on the proposed mechanism. The procedure of applying the static field is detailed in Section S2. Before presenting the results, it is helpful to review the behavior of electron drift velocity in graphene under an external electric field. In field-effect transistors, the drift velocity initially increases approximately linearly with the electric field, with a slope defined by the electron mobility. At higher electric fields, however, the drift velocity saturates, primarily due to the electron–phonon scattering [61]. To simplify our analysis, given this nonlinearity, we fix the drift velocity at v_x = 0.7v_F over a range of electric field strengths. The simulation results are presented in Fig. 3G to I. Under a relatively weak field E₀ = 1 × 10⁴ V/m, the increment of plasmons shows little deviation from the field-free case (Fig. 3G versus Fig. 3A). However, as the electric field strength increases to E₀ = 1 × 10⁵ V/m (Fig. 3H), a substantial increase in the growth rate is observed, accelerating the transition into the nonlinear stage. This confirms that a stronger static electric field markedly promotes the plasmon amplification. Furthermore, the corresponding plasmon spectra (Fig. 3I) reveal a gradual redshift in the fundamental mode frequency and an obvious enhancement in higher-order harmonics with increasing field strength.

Currently, these results are derived from calculations and analyses performed on ideal monochromatic perturbations. However, it remains to be verified whether the graphene plasmon can guarantee the purity of the spectrum during the THz amplification process in the presence of noise interferences. We define hybridization 1 as consisting of a monochromatic wave at 1.6 THz superimposed on random noise, and hybridization 2 as consisting of 2 monochromatic perturbations whose frequencies are 1.6 and 2.8 THz, respectively. In hybridization 2, the perturbation with 1.6 THz is considered as monochromatic noise. The signal-to-noise ratio (SNR) of hybridizations 1 and 2 are about 5 and 12.5 dB, respectively, according to the formula: SNR = 10lg (E_s/E_n), where E_s and E_n are the electric field amplitudes of the monochromatic perturbation and the noise perturbation, respectively.

From Fig. 4A to C, we can find that in the case of hybridization 1, the influence of noise on the spectrum during amplification is very limited, especially in the linear stage, where the effect due to noise is negligible. Figure 4D to F shows that the monochromatic noise in hybridization 2 can alter the spectrum markedly as higher harmonics are generated and amplified in the nonlinear stage. Especially when the electric field is saturated, the effect is remarkable, producing a superposition of 2 sets of THz frequency combs. Although monochromatic noise changes the results of the frequency comb in the nonlinear stage, it has little effect on the spectrum of the linear stage. From the results of the 2 cases, it is clear that the monochromatic noise with a higher SNR instead has a larger effect on the spectrum of the monochromatic perturbation than the random noise with a lower SNR. We argue that the random noise contains many frequency components. The components falling within the amplification region are limited and their amplitudes are not enough, so the effect on the spectrum of monochromatic perturbation is negligible. Compared to random noise, monochromatic noise contains only a single frequency component that falls within the amplification region, which makes it possible to alter the spectrum markedly. To summarize, we believe that the plasmon instability is highly resistant to interference and maintains the purity of the spectrum well before the electric field is amplified to saturation during THz amplification. This feature is advantageous for THz amplifiers.

We believe that the novel graphene plasmon instability has 2 potential applications. The first application is an ultra-broadband THz amplifier. By controlling the electron transit time to maintain operation in the linear stage, we can suppress the generation of higher harmonics, ensuring the amplification of the THz wave and the purity of its spectrum. For example, when we set the electron drift velocity v_x = 0.7v_F and the electron concentration n = 1 × 10¹² cm⁻² or n = 5 × 10¹² cm⁻², the linear stage time of graphene plasmon amplification is within 10 and 2 ps, respectively. Thus, the device’s operational length would be on the order of micrometers and submicrometers. The second application is the THz frequency comb generator, which requires a longer electron transit time to facilitate frequency multiplication and amplification in the nonlinear stage.

To experimentally verify the THz amplification mechanism investigated in this paper, here we analyze the challenges and the corresponding possible solutions in the following aspects.

Firstly, to achieve precise control over carrier concentration and drift velocity in graphene, the fabrication of graphene-based field-effect transistor structures with grating gates is a feasible option [37,46]. In these structures, the carrier concentration can be modulated by the gate voltage, and the source–drain voltage drives the electrons to accelerate until saturation. The output characteristic curve of the device characterizes carrier mobility and drift velocity. Secondly, realizing the proposed mechanism hinges on large scattering rate differences across the Fermi surface, which is easier to achieve at high drift velocities. Thus, fabricating graphene-based devices with ultrahigh carrier mobility is experimentally critical yet challenging. To achieve this objective, optimization can focus on 3 key aspects: monolayer graphene preparation, encapsulating dielectric selection, and electrode contact engineering. Among fabrication techniques for monolayer graphene, mechanical exfoliation and chemical vapor deposition (CVD) are dominant. The former offers high-quality samples at low cost, while the latter enables large-area synthesis [62]. Given the required micrometer-scale device dimensions, mechanical exfoliation is preferable for graphene preparation. For dielectric encapsulation, hexagonal boron nitride (hBN)/graphene/hBN heterostructures are ideal due to graphene’s ultrahigh mobility in hBN. This heterostructure exhibits low charge carrier inhomogeneity and ballistic transport over micrometer lengths at low temperatures [63]. Most studies fabricate such heterostructures via mechanical exfoliation and layer-by-layer transfer [64]. Critically, optimizing fabrication methods and electrical contacts has yielded devices with room-temperature mobility as high as 140,000 cm²/Vs [65]. This is achieved through van der Waals assembly of heterostructures, where 1D electrode contacts are formed by etching the top-layer hBN to expose graphene edges.

After fabricating the devices, exciting the graphene plasmon successfully is another challenge. In general, there are 2 approaches. The first is self-excited oscillation. Thermal noise distributed in the environment causes carriers in graphene to generate density perturbations, thus exciting the graphene plasmon [66]. In this way, 2 modes can be excited at the same time, i.e., downstream and upstream plasmons. The upstream plasmon can be amplified, while the downstream plasmon is attenuated when the plasmon instability studied in the manuscript is excited. The nature of this method lies in the fact that the thermal noise contains a large number of quasiparticles, which are plasmons with rich frequency information. Secondly, the graphene plasmon can be excited by external electromagnetic waves through momentum matching, which usually can be realized by gratings, nano-antennas, and near-field couplers [67]. This approach can be used to excite upstream plasmons by adjusting the coupling structure and the angle of incidence. It is also expected to produce coherent THz radiation.

To summarize, although it will be challenging to experimentally validate our results, some outstanding works prove the feasibility of carrying out the experiments to the best of our knowledge. This gives us strong expectations for future experimental observations of the mechanism by researchers.