Based on the above model settings, we simulate our device and evaluate the frequency spectrum of DC source-to-drain response voltage (dU). According to Dyakonov-Shur theory [20], dU is proportional to the intensity of THz signal, which, in turn, is proportional to the squared THz voltage amplitude. For a single-gate device, the response has the form [20]. d U = e V a m 2 4 m * s f ω (4) where ω is the angular driving frequency, f(ω) is a frequency-dependent function associated with the plasma wave (or damped electron wave) propagation properties. Generally, Eq. 4 applied to uniform-channel devices (i.e., single gate and no drain bias). Our recent work have shown that (4) can also be modified to evaluate the response in spatially nonuniform channel TeraFETs, such as exponentially-varying gate capacitance devices [2] and current-driven TeraFETs [47]. Besides, (4) does not take into account the device loading effect and thus only applies to devices with an infinite load impedance (e.g., an open-drain TeraFET). With a finite load resistance R _L, the output response voltage dU _o can be expressed by dU _o = dU _I/(1+ R _ch/R _L) [48, 49], where dU _I is the intrinsic response given by (4), and R _ch is the channel resistance. Due to the voltage divider effect, the loaded TeraFET has a lower response voltage compared to the unloaded ones under the same operating conditions.

The results under linear V _g(x) profile and 3 different N _s values under multi-gate structure are presented in Figure 2. These results are for a Si FET with 50% duty ratio (P _L = P _H = 50%, see Figure 1), and V _g0 = −0.2 V, V _th = 0.2 V, V _am = 2 mV, and at room temperature. Now the fundamental resonant frequency f ₀ is at 0.515 THz, and the device is driven into the subthreshold mode with a low electron density (∼10¹⁴ m^-2) but a relatively large voltage response [49, 50]. To improve convergence, the continuous-first-derivative transition regions are set between neighboring gate sections, and the relatively size of those regions (T _z, the ratio of total transition region size over the whole channel size) is fixed at T _z = 0.1 (see more details in Supplementary Material S1). The number of sections varied from 2 to 7. The voltage applied to different gates varied (see Figure 1B). Table 1 lists the major parameters used in the simulation.

Figure 2A shows the result under a linear varying gate voltage: V _g(x) = V _g0(1+α(x-0.5L)/L). We can see that with the increase of α (or the decrease of DC gate bias swing from source to drain since V _g0 is negative), dU decreases when f < f ₀, where f ₀ = S/4L is the fundamental resonant frequency, S is the plasma wave velocity [22, 50]. This region corresponds to the non-resonant operation region of the device. Using the methods in [2, 50], we can get the expression of DC response in this region (see Supplementary Material S1 for detailed derivations): d U = e V a m 2 4 m S 2 1 + β − 1 + β ⁡ cos 2 k r L cosh k 1 L cosh k 2 L (5) where β = 1/sqrt(1+(ωτ)⁻²), k ₁, k ₂ are wave vectors of the plasma wave: k 1 = a 1 * 2 S 0 S 2 + a 1 * 2 S 0 S 2 2 + i k 0 2 k 2 = a 1 * 2 S 0 S 2 + a 1 * 2 S 0 S 2 2 − i k 0 2 (6)

Here α 1 * = 1 V g 0 ∂ V g ∂ x , S 0 = e V g 0 m * , S = η e V t m * 1 + exp − U 0 η V t ln 1 + exp U 0 η V t . Besides, k ₀ = k ₀ = (ω/S ² τ)^0.5 is the wave vector in the uniform channel (α ₁ ^* = 0), k _r is the real part of k ₁ or k ₂. A transition of variation trend with respect to α occurs at around f = f ₀. Beyond f ₀, the plasmonic resonance can be achieved, and dU decreases with increasing α. Now the response curve does not follow Eq. 6. Within α ∈ [0,0.5], the maximum improvement of dU is around 20%. Those results agree with our observations of linearly varying gate capacitance or threshold voltage [2].

Figure 2B shows the result of dU vs α under N _s = 3. Compared to Figure 2A, the 3-segment multi-gate TeraFET exhibits a distinct response profile, and Eq. 6 fails in this case. As f rises, the response voltage oscillates, and the variation trend of dU with respect to α changes multiple times. If we define the regions where dU increases with rising α as the “enhanced” regions, and the regions where dU decreases with rising α as the “suppressed” regions, we can see that the enhanced and the suppressed regions appear alternatively with the increase of frequency. More interestingly, the positions of those regions are directly related to the number of gate sections. For example, the peak response voltage in the first enhanced region (which is also the peak dU in the whole frequency range) is at f = 3f ₀, the position of response valley in the first suppressed region is f = 6f ₀, and the position of peak response in the second enhanced region is at f = 9f ₀. Thus, we conclude that the frequency at which the maximum response is reached is around f p = N s f 0 (7) and the frequency gap between two adjacent peaks or valleys is d f p = 2 N s f 0 (8)

Eqs 8, 9 indicate that we can selectively enhance the resonant detection responsivity in a given frequency band by properly selecting the number and length of gate sections in a TeraFET. This result is similar to those reported in [34, 51], where the THz transmission spectra was controlled by the gate separations in grating-gating graphene FETs, and the resonant frequency was determined by the unit finger gate width (∼L/N _s). However, unlike the current-driven boundary condition used in [34, 51], here we used an open-drain condition. Thus the enhancement and suppression of voltage response are independent of plasmonic instability or DC field effects, and solely originate from the setup of plasmonic cavities. Besides, the TeraFETs used in this work operate in the subthreshold region, in which the electron density and plasma wave velocity are less sensitive to the gate bias variations [50]. This helps stabilize the resonant frequency at each order (since f ₀ = S/4L, where f ₀ is the fundamental resonant frequency). This explains why there is no significant redshift or blueshift observed in Figure 2 as α changes.

Eqs 8, 9 can be further verified by the simulations under other N _s values. For example, in Figure 2C where N _s = 4, the peak frequency is at 4f ₀ and the distance between two adjacent peaks or two adjacent valleys are 8f ₀. In Figure 2D where N _s = 5, the values of f _p and df _p are 5f ₀ and 10f ₀, respectively. Also, a 100% increase in dU (compared to the uniform channel case) is achieved when α reaches 0.5. Note that the peaks and valleys are not located at the fundamental resonant frequency, but at the higher order harmonics. Therefore, the introduction of multiple gate sections activated the harmonic components in the system, resulting in the distribution of enhanced and suppressed regions. The underlying mechanism could be related to the reflection of plasma waves or carrier drift between neighboring sections due to the carrier concentration barriers. Those reflections change the wave propagation properties (i.e., k ₁ and k ₂) and shorten the effective channel length, thereby leading to the excitation of harmonic peaks and valleys. Figure 3A shows the spatial distribution of gate-induced field (dV _g/dx) along the channel. The abrupt change of DC gate bias in the narrow transition regions creates a large field on the order of 0.1 MV/cm. The electrons passing the transition regions get accelerated or de-accelerated, forming the separated velocity distribution regions, as demonstrated by the velocity distribution contour plot in Figure 3B. Each curve in Figure 3B represents the carrier velocity distribution u(x) at a given moment in one AC period, and 50 consecutive moments are included. The separated velocity distribution regions could induce the reflections of plasma waves in between, thereby altering the DC response properties.

The above harmonic excitation mechanism can be seen as a result of abrupt changes in channel properties, as opposed to the gradual changes reported in our previous work [2]. In a gradually varying channel, the response performance is related to the changing rate of channel parameters (e.g., the gate capacitance, threshold voltage, DC gate bias). While in multi-gate setup, we can verify from simulation that the response dU is insensitive to the transition region size T _z (see supplementary material). This indicates that the response profile is now level-sensitive, as opposed to the gradient-sensitive ones in [2]. Therefore, the analytical approaches developed in [2] can no longer be applied here.

To further investigate the variation trend of dU with frequency, we define a differential response voltage ddU = dU(α = 0.3) - dU(α = 0), and plot its frequency profile at different Ns values, as shown in Figure 4. Here ddU signifies the net enhancement or suppression of dU at α = 0.3 as compared to the uniform channel case. In Figure 4A, the amplitude of ddU rises with the increase of N _s. This suggests that the enhancement effect strengthens as the channel becomes more non-uniform. With the rise of frequency, ddU oscillates and exhibits multiple peaks and valleys, as shown in the separated plots Figures 4B–D. For quantitative analysis, we plot Figure 5 where ddU _max, f _p and df _p as functions of N _s are presented. One can check that the f _p and df _p curves follow Eqs 8, 9. The ddU _max increases with the rise of N _s, but a saturation trend is observed when N _s becomes large. This saturation could be related to the change of wave reflection characteristics as the length of each gate section shortens, which sets a limit to the maximum improvement of dU.

To better understand how the response changes with frequency and gate structure, we create a map of ddU in a N _s-f/f ₀ scale, as shown in Figure 6. In the map, the enhanced regions are exhibited as “mountains” while the suppressed regions are presented as “valleys”—a result of the present ddU definition. The highest mountain group is located at f = N _s f ₀, as shown in Figure 6A, which corresponds to the maximum (the first) resonant peak in each case. The second mountain series are at f = 3N _s f ₀, demonstrating the secondary resonant peaks. Between these two mountain groups is a valley group located at f = 2N _s f ₀. In general, the mountain clusters can be expressed by f = (2n+1)N _s f ₀, where n = 0,1,2 … , and the valley clusters follow f = (2n+2)N_sf₀ .

Figure 6B shows the direct comparison of the heights of different mountains (i.e., the amplitudes of response peaks). Clearly, the mountain height in each group increases with the increase of N _s, and the average/maximum height in the first mountain group is much larger than that in the second mountain group. Thus, to achieve a high response, the TeraFET should operate in the first mountain group, and in general a large gate section number is preferred.

The results in Section 3.1 and Section 3.2 demonstrate that adopting periodic multi-gate structure in TeraFETs can effectively alter the DC voltage response and achieve over ∼100% improvement in dU at certain frequencies. The amplitude of dU can be tuned by N _s and α. In general, a larger N _s or α leads to a higher responsivity in the enhanced region, but the values of N _s or α cannot grow infinitely due to several built-in limits. Here we discuss those limits.

1) Breakdown voltage (vertical). To prevent the breakdown of the barrier material, the following is required

For example, if E _b = 3 V/nm, d _b = 4 nm, |V _g0| = 0.2 V, we get α < 46.

2) Breakdown voltage (transverse). Let D denotes the transition region length between two gate sections, and D is related to T _z. To prevent dielectric breakdown in the transition region, we need

If E _b = 3 V/nm, |V _g0| = 0.2 V, D = 2 nm, we get α < 30.

3) Conductivity limit. When the gate bias decreases in the subthreshold region, the carrier concentration can reduce to very low, so as to choke the current conduction. Assume that the minimum conductivity required for sustaining current conduction is σ _cr = eμn _cr, where n _cr is the critical carrier density. Using Eq. 3, we get:

If n _cr = 10¹³ m^-2, V _g0 = −0.2 V, V _t = 0.026 V (T = 300 K), η = 4, we get |α| < 1.29.

4) Process limit. The fabrication lab conditions determine the maximum number of separated gates that can be built in a TeraFET. If the minimum achievable size is L _min, then we get N _s-max = [L/L _min], where [k] denotes the nearest integer that does not exceed k. For example, with L = 250 nm, L _min = 65 nm, we get N _s-max = 3.

The above conditions, along with other more delicate mechanisms (e.g., the self-capacitance and the built-in voltage between two adjacent gate sections), set limit to the tuning of dU in periodic multi-gate TeraFETs. Despite all those constraints, an ∼100% improvement can still be achieved near the maximum resonant peak, as demonstrated in Figure 2D.

Table 2 summarizes the THz detection performance of TeraFETs with non-uniform gate capacitances and periodic multi-gate structures. The data in this table are taken either from [2] or from this paper. It can be seen that using gate capacitance profiling, the maximum reported response improvement was around 20%. With periodic multi-gate design, the response improvement can rise beyond 100% if sufficient gate sections are introduced. This further exhibits the advantage of periodic multi-gate design over gate capacitance profiling.