Modeling of Antenna-Coupled Si MOSFETs in the Terahertz Frequency Range

We report on the modeling and experimental characterization of Si complementary metal-oxide-silicon (CMOS) detectors of terahertz radiation based on antenna-coupled field-effect transistors (TeraFETs). The detectors are manufactured using Taiwan semiconductor manufacturing company (TSMC's) 65-nm technology. We apply two models—the TSMC RF foundry model and our own advanced design system (ADS)-hydrodynamic transport model (HDM)—to simulate the Si CMOS TeraFET performance and compare their predictions with respective experimental data. Both models are implemented in the commercial circuit simulation software keysight ADS. We find that the compact model TSMC RF is capable to predict the detector responsivity and its dependence on frequency and gate voltage with good accuracy up to the highest frequency of 1.2 THz covered in this study. This frequency is well beyond the tool's intended operation range for 5G communications and 110-GHz millimeter wave applications. We demonstrate that our self-developed physics-based ADS-HDM tool, which relies on an extended 1-D HDM and can be adapted readily to other material technologies, has high predictive qualities comparable to those of the foundry model. We use the ADS-HDM to discuss the contribution of diffusive and plasmonic effects to the THz response of Si CMOS TeraFETs, finding that these effects, while becoming more significant with rising frequency, are never dominant. Finally, we estimate that the electrical noise-equivalent power (perfect power coupling conditions) is on the order of 5 pW/$\sqrt{\mathrm{Hz}}$ at room-temperature.


Introduction
Antenna-coupled field-effect transistors (TeraFETs) have emerged as powerful devices for a variety of terahertz (THz) applications such as frequency multiplication and mixing [1,2,3,4,5,6], as well as direct THz power detection [1,7,8,9,10,11,12,13,14,15,16,17,18,19].The detectors can be readily integrated in electronic circuits and provide sufficient sensitivity at room temperature which makes them well suited for applications.Among these, TeraFETs fabricated in Si CMOS technology have reached impressive performance levels [13,20] with a high sensitivity over almost the entire THz frequency range.With the growing variety of possible applications of field-effect transistors (FETs) in the THz range, the need for reliable design and modeling platforms for such devices has increased significantly.Here, we focus on detector simulation and experimentally validate our results.In [21] we presented our simulation approach for TeraFETs -acronym: ADS-HDM.The ADS-HDM is implemented into the circuit simulation environment Keysight ADS -an industry standard.Its implementation is based on an extended version of the one-dimensional (1D) hydrodynamic transport model (HDM), which has been used in [1].The frequency-dependent channel impedance is determined from the self-consistent solution of the HDM equations.This physics-based model accounts for the mixing process in the FET channel involving the charge-carrier density waves (plasma waves).They are excited by the incident THz radiation via an antenna element which couples the radiation asymmetrically to the FET channel, either via the source-gate or the drain-gate port [22,13,23,24,25].The evidence for the existence of plasma waves (which are overdamped in the case of Si CMOS at room temperature) has been strengthened by experimental findings described in [23,26,27,28].In this work, we present and discuss two simulation approaches to simulate the detector performance of Si CMOS TeraFETs, our in-house tool ADS-HDM as well as the TSMC RF foundry model (TSMC RF) of the Taiwan semiconductor manufacturing company (TSMC).Prior to the simulations we propose an approach to determine the exact THz input power present at the transistor terminals when the detector is illuminated with a Gaussian intensity profile.The model predictions are compared quantitatively with the optical voltage responsivity determined between 0.4 and 1.2 THz and put into perspective with simulation results obtained with the TSMC RF.The ADS-HDM tool allows to switch the plasmonic effects on and off.Using this feature, we quantify the plasmonic contribution to the THz detector responsivity.Finally, we discuss the upper sensitivity limits of the 65-nm Si CMOS TeraFETs.

Hydrodynamic Model description
In order to simulate the excitation of density waves of the charge carriers in the channel of the Si MOSFET detectors, we use the following 1D hydrodynamic transport equations consisting of the continuity equation Eq. 1a and the momentum balance equation Eq. 1b, both written in terms of the current density j = −qnv, where q is the electron charge, n = n(x) represents the carrier sheet density in the FET's inversion layer along the channel (x-direction), and v = v(x) stands for the local mean carrier velocity.Note thatcontrary to the model presented in [1,22] -, the second term on the right hand side of Eq. 1b also takes carrier diffusion into account.Other quantities of the equations are the longitudinal electric field E x = −∂ x φ as determined by the channel potential φ(x), the electron effective mass (m * = 0.26 m 0 for electrons in Si in units of the free-electron mass m 0 ), and the momentum relaxation time τ p of the electrons.The temperature of the crystal lattice is assumed to be room temperature (T L = 293 K).k B denotes the Boltzmann constant.The dependence of n(x) on the externally applied gate voltage is taken into account by the unified charge control model (UCCM) for n-channel MOSFETs [29] n where C ox = ϵ R ϵ0 d is the channel capacitance per unit area for a thickness d of the gate insulator.For d, we assume a value of 3.05 nm. 3 ϵ R = 3.9 is the relative permittivity of the gate insulator, and ϵ 0 the vacuum permittivity.U T and U th = k B T L q are the threshold voltage and the thermal voltage, respectively.γ represents a non-ideality factor, which determines the subthreshold slope of the drain current with regards to the gate-to-source voltage U GS .The gate-to-channel voltage U gc (x) = U GS − φ(x) is given by the difference of U GS and the channel potential φ(x), the latter being determined self-consistently by jointly solving Eqs.1a, 1b and 2. The numerical evaluation of this set of equations is performed in the circuit simulation environment of Keysight ADS using a novel distributed-element method approach to integrate the transport equations as presented in [21] for AlGaN/GaN HEMTs.The steady-state solution is obtained by the harmonic balance technique.A schematic of the ADS-HDM implementation is shown in Fig. 1(a).A significant advantage of this simulation approach versus recently developed TCAD models for the THz frequency range [31,32,33,34] is that in the case of circuit based models like ADS-HDM [21] or recently developed compact SPICE models [35,36], extrinsic and parasitic device components are readily included in the full simulation domain (as indicated in Fig. 1(a)), while for TCAD simulations they have to be included via more complex mixed mode simulation methods [37].At the same time, the full physical model for the description of the plasma-wave mixing process in the transistor's channel is retained at the model's core (hydrodynamic domain).By manipulating or switching off individual terms in Eqs.1a, 1b and 2, one can study and identify the influences of the various physical processes included in the equations.A number of input values for the parameters of the simulations, as presented in Table 1, can be obtained from the measured DC drain-to-source resistance (see below).This makes it easier to test whether variations of the THz performance of each TeraFET may be caused by fabrication-related issues which also affect the DC device properties.The blue box (full simulation domain) represents the THz detector model (is taken in the same manner for the TSMC RF foundry model), here for source radiation coupling conditions.At its core (hydrodynamic domain), the hydrodynamic transport equations are used to calculate the rectified signal of the detector on the basis of plasma-wave-modified resistive mixing in the transistor's channel.The equations are solved in the Keysight ADS simulator itself.The FET parasitics are taken into account in the FET circuit domain.Parasitic components considered are the source and drain overlap C ov and fringe capacitances C f r , the contact resistances R C,S and R C,D , and the gate resistance R G .The antenna parameters including the antenna impedance Z ant are estimated from EM wave simulations (CST Studio Suite).For details of the HDM implementation, we refer to [21].(b) Cross-sectional view of the detector implementation in TSMC 65-nm process.The patch antenna is implemented underneath a passivation layer in the top metal layer MT.The ground plane is placed in lowest to metal layers M1 and M2.

65-nm Si CMOS TeraFETs
The detectors studied in the following are designed for a commercial 65-nm complementary metal-oxide-silicon (CMOS) foundry process (CLN65LP, type 1.2V_mLow_Vt_MOS) of the Taiwan Semiconductor Manufacturing Company (TSMC).Fig. 1(b) displays a cross-sectional schematic view of the detector design.In brief, a narrow-band rectangular patch antenna, which is placed in the last metal layer M10, collects the THz radiation and guides the signal through a via at one of the narrow sides of the patch to the source of a single transistor placed directly under the edge of the patch.The height of the patch above the ground plane (metal layer M1 and M2) is 8.41 µm.In the following, four different devices are studied with respect to their detector performance.The devices differ only in terms of their respective patch length L patch (see micrograph with the values in Fig. 2(a)), which results in different resonant frequencies.The patch width is fixed to 40 µm.A 1.3-µm-thick passivation layer covers the devices.An important figure of merit, which determines the efficiency of our detectors, is the optical voltage responsivity R op V . 4It is defined as the rectified voltage (or current) signal between the drain and source terminal ∆U DS with regards to the  9) numerically (solid lines).In both cases the FWHM of the Gaussian beam is determined from Eq. ( 7) (as shown as black line in (e)).
available optical THz input power P op in , the total power of the incident THz beam in free-space [38,39]: In this work, the responsivity is measured with the help of an optoelectronic spectroscopy system (TeraScan 1550, Toptica Photonics AG, Graefelfing (Munich)), where the optoelectronic detector is replaced by our TeraFETs.∆U DS is the beam size.It is calculated as the optical responsivity multiplied by the ratio of the (measured or calculated) beam cross-sectional area at the antenna to the (calculated) antenna cross-sectional area [38].
measured using the lock-in technique.For this purpose, the bias to the photoconductive emitter is modulated with a sine wave at f mod = 7.62 kHz.A schematic drawing of the measurement setup is displayed in Fig. 6(a).Continuous-wave THz radiation between 0.1 and 1.2 THz is generated by an InGaAs photomixer emitter.Its radiation is guided to the Si CMOS TeraFET via four off-axis paraboloid (OAP) mirrors with 90°-beam-deflection geometry.The THz beam path is purged by dry nitrogen in order to minimize radiation absorption by water vapor.The total power P op in of the incoming THz beam in front of the Si CMOS TeraFET, placed in the focal point of the last OAP mirror, is measured using a calibrated large-area Golay cell.
The peak-to-peak rectified voltage between drain and source contacts is obtained as ∆U DS = α LIA • U LIA , where U LIA is the lock-in amplifier's output signal.It is proportional to the RMS of the first Fourier fundamental component of the detector signal.α LIA represents a lock-in correction factor.For an ideal sinusoidal modulation (e.g by electronic chopping) of an emitter, the lock-in correction factor is found to be a constant value of α LIA,corr = 2 √ 2. However, inspection with an oscilloscope shows that the THz wave is not modulated sinusoidally, but rather exhibits a more square-wave-like waveform with a detailed shape which depends on the beat-note frequency of the laser radiation.Measured results and the derivation of the THz-frequency-dependent values of α LIA,corr are discussed in Appendix A. α LIA,corr is found to increase linearly with frequency from ∼ 2.4 (at 0.1 THz) to ∼ 3.1 (at 1.2 THz).

Power coupling theory
For a proper comparison of the simulated optical voltage responsivity with the measured one to be addressed later, the electrical THz input power P el in,t at the source transistor terminal (here for source coupling conditions) has to be estimated.For this purpose, we need to consider multiple losses suffered by the incoming radiation before it is converted into the electrical input signal.P op in denotes the power of the incoming THz beam directly before the passivation layer.The electric power P el in,ant , which one obtains from this optical power at the antenna is calculated in the following way: Here, η op takes into account all possible optical losses of the incident THz radiation on the path from the source to the detector, which are not considered by the electromagnetic (EM) wave simulations.For example absorption and reflection losses associated with the substrate and a lens, which may be attached to it [25,39].In this work, we set η op = 1, since the incident THz light is coupled directly from free space to the antenna as indicated in Fig. 1(b), neither passing through a substrate nor a substrate/superstrate lens (backside illumination is obstructed by the ground plane).η scat represents the absorption efficiency of the antenna.It includes the residual scattering losses that occur in a receiving antenna element even in case of a perfectly matched load.For large planar antennas it is known to be η scat ≈ 0.5 [40,38].The final loss factor is the frequency-dependent Gaussian beam coupling efficiency η gauss (ν).From a geometrical point of view, η gauss (ν) considers collection losses resulting from the discrepancy between the spatial extension of the incident Gaussian beam's intensity profile, given by [41] and the effective area of the antenna [42] A Here z R = πν 4c0 • D f 2 is the Rayleigh range and k 0 = 2π λ0 = 2πν c0 the wave vector in free-space with speed of light in vacuum c 0 .I 0 denotes the intensity in the center of the Gaussian beam.D f represents the minimum 1/e 2 Gaussian beam diameter where d L is the diameter of the collimated optical beam before the focusing element and f L is the focal length of the focusing element (in this work an OAP mirror with f L = 2" = 50.8mm).The total incident THz power is distributed over I gauss (r, φ) such that A ef f of the receiving antenna element is determined in the direction of maximum gain G. 5 ϵ ant is the antenna efficiency and D the antenna directivity.Both, A ef f and ϵ ant , can be obtained from EM simulations (CST Studio Suite, shown in Fig. 2(d)) or from direct measurements (e.g.absolute gain methods or gain transfer methods).For more details about the different methods to estimate the effective area of a THz detector, we refer to [38].Our EM wave simulations take the full three-dimensional structure of the Si CMOS TeraFET into account (i.e., all structures shown in the cross-sectional view of Fig. 1(b)).This includes the vias from the antenna to the transistor and the ground-plane, the surrounding metallization, the different dielectric layers and the passivation layer.The lumped element port is placed at the position of the transistor.In general, the EM wave simulations of ϵ ant (ν) and thus also of η gauss (ν) consider both conduction losses of the metallization and possible dielectric losses due to the surrounding dielectric layers.Consistent with the real device, the metallization is simulated as lossy copper, so conduction losses are considered in our simulations.Dielectric losses are not considered here because the dielectric layers are made of insulator materials with zero loss tangent.We want to clarify that the Gaussian beam coupling efficiency η gauss (ν) is closely related but not identical to the gaussicity [43] -a quantity typically employed for substrate-lens-assisted antennas.These typically are configurations, where a dielectric lens (e.g.hyper-hemispherical silicon lens) is placed onto the detector's substrate to enhance the Gaussian beam coupling efficiency under backside illumination conditions [44,45,43].The gaussicity can be estimated from the overlap integral of the normalized simulated antenna pattern and the best focused Gaussian beam in the angular domain [43,39].This calculation imposes that the effective area of the antenna is larger than (or of the same size as) the Gaussian beam diameter at the position of the antenna plane (after the dielectric lens and the substrate).In this work, we determine η gauss (ν) for free-space optical coupling conditions, where the latter condition is not fulfilled.In the following we assume that the detector is positioned in the focal point of a focusing element such as an OAP mirror (see Fig. 6(a) in Appendix A) and that the evolution of the optical beam in the experiment can be treated on basis of Gaussian beam theory [41].Under these circumstances η gauss (ν) can be calculated for arbitrary frequencies from where FWHM = ln(2)/2 • D f is the frequency-dependent minimum full-width at half maximum of the incident Gaussian beam given in the antenna plane (in focal point of the OAP mirror) and A spot = π • (FWHM/2) 2 is the spot size of the beam.Note that the integral for β(ν) has no analytical solution and needs to be solved numerically.However, we find that β can be approximated by β ≈ π/4 − ln(2) • A eff /A spot , which leads to η gauss ≈ ln(2) • A eff /A spot .For details on the derivation of the equation set Eq. ( 9), we refer to Appendix B. Precise knowledge of d L is therefore required to determine the Gaussian beam intensity profile at the detector position.
For this purpose, we estimate the FWHM (or D f ) at the focal point of the OAP mirror by direct measurements of the Gaussian beam intensity profile I gauss (x, y) using an automated X-Y micrometer translation stage as shown in Fig. 2(b) at the respective resonant frequency of each detector A1-A4.We want to emphasize that I gauss (x, y) may have to be de-embedded when using highly directional antennas with a narrow beamwidth.The 3-db beamwidth (angular width) of the patch antennas in this work is roughly 90°over a wide range of frequencies.From I gauss (x, y) we determine the FWHM of the Gaussian distribution alongside X (shown as red line) and Y (shown as green line) using a Gaussian fit model (see Fig. 2(e)).Then the mean FWHM can be obtained from FWHM = FWHM x • FWHM y (shown as green line).Consequently, we take the observed frequency roll-off of the measured FWHM to define d L in our measurement system.Best agreement with direct measurements of the FWHM is observed for d L = 30 mm (calculated from Eq. ( 7), shown as black solid line in Fig. 2(e)).
In Fig. 2(f) we plot the calculated Gaussian beam coupling efficiency η gauss (ν) obtained from the numerical solution of Eq. ( 9) for each patch antenna (A1-A4) taking the simulated effective areas A ef f and the experimentally determined minimum FWHM into account.Besides that, we compare the solution of Eq. ( 9) with values of η gauss obtained from the direct evaluation of the numerical area integral for P ant (see Eq. ( 13) in Appendix A) and P gauss (Eq.( 8)) using (i) the measured and (ii) an ideal symmetric Gaussian beam intensity profile I gauss (x, y).The latter is determined from the minimum FWHM (Eq.( 7)) at the respective frequencies.From Fig. 2(f), even for the A3 and A4 with larger effective antenna area, only about 10% of P op in can be coupled into the antenna element in the best case.Finally, the where η m (ν) accounts for the impedance mismatch of the antenna impedance Z ant (ν) (see Fig. 2(c)) to the transistor channel impedance Z t (ν).Best electrical power coupling to the transistor is ensured whenever conjugate impedance matching conditions are present Z ant = Z * t .Note, that impedance mismatch effects (η m (ν)) are inherently embedded in the circuit simulation software ADS (e.g. by taking a Power Source at a certain frequency with complex impedance) and therefore, for the simulations presented here, we use P el in,ant as THz power input parameter.Besides the simulated antenna parameters Z ant (ν), A ef f (ν) as well as the Gaussian beam coupling efficiency η gauss (ν) all other unknown input simulation parameters, in case of ADS-HDM, are obtained from DC drain-to-source resistance (R DS ) measurements.The measured resistance curves are fitted to a simple serial resistance model for the MOSFET channel, which is based on the UCCM.The extracted simulation parameters are shown in Table 1. 6

Simulations and Comparison with Measurements
In order to validate the implementation technique and the fitting routine discussed previously we compared the simulated R DS as a function of U GS with measurements.In Fig. 3, R DS is shown for TSMC RF and ADS-HDM together with the respective measured data for a fixed drain-to-source voltage of U DS = 10 mV.We observe good agreement between the ADS-HDM and data, which is a direct consequence of the specific simulation parameter extraction performed for each device separately (see Table 1).Contrary TSMC RF exhibits only near-quantitative agreement with data close to U T .In the inversion regime (U GS >U T ) the measured R DS are slightly higher than predicted by TSMC RF.To the best of our knowledge, the discrepancy originates from globalized simulation parameters of the TSMC RF foundry model.These globalized simulation parameters are only adapted depending on the CMOS fabrication process (here CLN65LP, type 1.2V_mLow_Vt_MOS) and the device dimensions (length and width of the channel).Therefore, the foundry model can not account for possible variations (e.g.due to imperfections, defects etc.) of the electrostatic properties within a foundry fabrication run.A comparison of the frequency-dependent measured optical voltage responsivity (top panels) of devices A1-A4 with the respective simulations -TSMC RF and ADS-HDM -can be seen in Fig. 4. The final unknown simulation parameter of the ADS-HDM -the total parasitic capacitance per micron channel width C ov+f r,W −1 , which includes the overlap C ov and fringe C f r capacitances -is determined from the comparison of the measured responsivity roll-off of the detector characterized at the highest THz frequencies (A4) with the respective ADS-HDM simulations.Best agreement is found for C ov+f r,W −1 = 0.75 fF/µm.Similar values have been reported for Si CMOS TeraFETs manufactured with the 150-nm technology by TSMC [13].Besides that, in simulations of the detector voltage response ∆U DS , the calculation of the DC read-out is performed against an effective loading impedance [46], taking into account the lock-in amplifier's input impedance R in (10 MΩ) and the effective capacitance of the cable (C L = 0.5 m • 0.95 pF/m).
We obtain a near-quantitative agreement between both models and measurements at low THz frequencies (Device A1 and A2).Both models predict magnitude, peak position and the roll-off of the measured signal quantitatively.Only far away from the peak (at 0.4 THz for A1 and A2) we observe a discrepancy between data and simulations.As presented in Fig. 2(c) the EM wave simulations of the antenna impedances Z ant suggest an additional antenna resonance at  -A4).U DS is set to 10 mV for all measurements and simulations.The dashed vertical line represents the extracted U T for the respective devices using a fit model based on the UCCM (see Table 1).For both models we assume a fixed residual resistance of R p,ext = 300 Ω for a short-circuit protection FET, which is connected in series with the drain electrode of all respective transistors.∼ 0.35 THz, which is not observed in the THz characterization measurements.At higher THz frequencies (Device A3 and A4) the predicted peak position and magnitude of the optical voltage responsivity of both models starts to differ.The observed divergence can be related mainly to different predicted losses due to imperfect conjugate impedance matching (Z ant ̸ = Z * t ) between Z ant and Z t .The latter is determined by either TSMC RF or, in case of ADS-HDM, by the hydrodynamic transport equations (Eq.( 1)) together with the UCCM (Eq.( 2)).In the bottom panels of Fig. 4 the real part and imaginary part of the simulated complex conjugate transistor channel impedance Z t * is depicted together with Z ant for each device.In addition we show the calculated impedance mismatch factor η m (Eq.( 10), shown as dots, right axis).Both models suggest that -especially in case of device A3 -perfect conjugate impedance matching (η m ∼ 1) is reached.Remarkably both models predict very similar Z t .Despite the fact that TSMC RF has not being adapted specifically for the THz frequency range, it exhibits good agreement with measurements and with ADS-HDM.In addition to the frequency-dependent measurements we also measured and simulated the optical voltage responsivity R op V as a function of gate voltage U GS .In Fig. 5(a), we compare the measured optical voltage responsivity of devices A1 to A4 with the results of simulations with ADS-HDM and the TSMC RF for fixed excitation frequencies.Again, we observe good agreement between both models and measurements.Best agreement is obtained in the inversion regime above U T (marked as vertical dashed line) of the respective transistor, while agreement gets worse below U T , where R op V (U GS ) peaks.The deviations are attributed to unknown details of power-coupling and U GS -dependent impedance matching.In addition to results for the full-fledged ADS-HDM, we also plot the predicted purely resistive-mixing response of the transistor, where the diffusive and plasmonic terms of the HDM are neglected and only the low-frequency limit is considered (∂ t j ≪ j/τ p ), such that Eq. (1b) is reduced to j = σ • E x , with the conductivity σ = (q 2 τ p n)/m * .The difference between the full ADS-HDM and the purely resistive self-mixing model increases with frequency.The responsivity as predicted by the full model is larger than that of the purely resistive mixing, albeit not by a fixed factor for all values of U GS .The maximal difference reaches values of about 20% for device A4 at 1070 GHz and U GS = 0.4 V.This increase is evidence that diffusive and plasmonic features cannot be neglected as the frequency rises.Their influence, however, does not become dominant up to the highest frequencies of our studies, and an approximative treatment of rectification by purely resistive mixing yields fair results.The presented results of measurements and simulations of the optical voltage responsivity have only considered the measured detector response with respect to the available optical THz input power P op in .For a detector sensitivity, noise contributions to the detected signal must be taken into account.The most important figure of merit in this context is the Comparison of the simulated frequency-dependent real (solid lines) and imaginary (dashed lines) antenna impedance Z ant (green) with corresponding simulations of the complex conjugate transistor channel impedance Z t * using TSMC RF (blue) and ADS-HDM (red) for the different Si CMOS TeraFET devices.Estimated impedance mismatch factor η m (dots, right axis) is calculated from Eq. ( 10) for TSMC RF (blue) and ADS-HDM (red).noise-equivalent power (NEP).The NEP is defined as the ratio of the detectors voltage or current noise spectral density (V N or I N ) in a 1 Hz bandwidth (∆f ) and the detectors voltage or current responsivity, respectively.It can be defined with regards to the optical or electrical detector responsivity (R op V or R el V ).For measurements of the voltage response we have where the voltage noise spectral density V N , in case of zero drain-to-source bias (U DS = 0 V) during THz detection operation, is dominated by Johnson-Nyquist (thermal) noise, such that [47,48].This assumption is based on noise characterization measurements of Si CMOS TeraFETs manufactured by TSMC's 90-nm technology [48].In Fig. 5(b) we present a comparison of the measured optical NEP with simulations obtained by the TSMC RF and the ADS-HDM.Both models correctly predict the overall magnitude and gate voltage dependency of the measured NEP op V close and above the U T .Below threshold and at higher THz excitation frequencies (A3 and A4) the ADS-HDM achieves better comparability with data.In order to show possible improvements of the detector  V of the respective devices.TSMC RF (blue filled squares) and ADS-HDM (red filled squares) simulations of the voltage-dependent electrical NEP (NEP el V ) with regards to the maximum available power delivered from the antenna P el in,ant .In addition, the simulated electrical NEP calculated with regards to the available power at the source transistor terminal P el in,t is depicted for TSMC RF (blue dashed line) and ADS-HDM (red dashed line).
performance (e.g. by enhancing the Gaussian beam coupling efficiency η gauss or reducing the back-scattering effect of the antenna) we present simulations of the electrical NEP el V , which represents the intrinsic lower limit of the detector sensitivity.We perform this calculation for two cases, where the electrical responsivity is referred (i) to the available power originating from the antenna element R el V = ∆U DS /P el in,ant and (ii) to available power at the Source transistor terminal R el V = ∆U DS /P el in,t .In the latter impedance mismatch losses η m are excluded.These losses play a more dominant role at lower THz frequencies, where the device impedance Z t changes more strongly with applied gate voltage.As already indicated by the estimated values of η gauss (see Fig. 2(f)) and also visible in Fig. 5(b) there is more room for improvements of the detector sensitivity for A1 and A2 -with respective small effective antennas areas (see Fig. 2(d)) -then for A3 and A4.For example, if η gauss is increased to ∼ 1 by reducing the FWHM of the Gaussian beam intensity profile in the antenna plane using common methods such as dielectric lenses (placed here on the patch antenna), one can expect that the detector performance increases by a factor of 40 for A1 and 12 for A4.In Table 2 the best determined optical NEPs for the respective detectors A1-A4 are shown together with the predictions by the TSMC RF and the ADS-HDM at the respective frequencies and gate voltages.In addition we present the simulated electrical NEP (with regards to the P el in,ant ) of the respective devices at the same frequency and gate voltage conditions.The simulations suggest that a minimum NEP of 5 − 6 pW/

√
Hz could be achieved at room temperature in ideal optical

√
Hz if the backscattering effect of the antenna is taken into account (η op = 1, η gauss = 1, η scat = 0.5).This in close agreement to the purely analytical calculations of the electrical NEP for 150-nm Si CMOS TeraFETs as presented in [13].

Conclusion
In summary, the detector voltage responsivity and sensitivity of four 65-nm Si CMOS TeraFETs in the THz frequency range of 0.4 and 1.2 THz and the gate voltage range of 0 to 1 V were investigated in detail using two transistor transport models -the TSMC RF foundry model and our in-house ADS-HDM circuit implementation.The THz characterization measurements were performed under nitrogen atmosphere and the accurate peak-to-peak response of our detectors was determined from time-domain data acquisitions in order to ensure a suitable database for both models.The amount of THz power injected into the transistor terminals by the Gaussian beam via the antenna was calculated based on geometric constraints.A near-quantitative agreement between model and experiment was observed over a wide range of THz excitation frequencies and gate voltages.The simulations indicate that the 65-nm Si CMOS TeraFETs perform close to perfect conjugate impedance matching and further improvement of their detector performance down to ∼ 10 pW/

√
Hz at room temperature can be expected if the efficiency of Gaussian beam coupling is improved.By comparing the two models with each other and with experimental I/V and THz characterization, we find that the TSMC RF foundry model is applicable beyond its actual design frequency range for 5G communications and 110 GHz millimeter-wave applications.It can be used with confidence to make predictions about detector performance in the investigated operating range of our Si CMOS TeraFETs up to 1.2 THz.However, the experimental data of A3 and A4 and their comparison with the simulations suggest that the TSMC RF model may lose its validity at frequencies above 1.2 THz at which the plasmonic mixing terms in Eq. (1b) are likely to begin to dominate the device response.Our work extends recent advances in compact SPICE circuit modeling of Si MOSFETs in the THz frequency range [35] by two more simulation approaches.In this work, we compare our simulations directly with experimental data and achieve near-quantitative agreement over a wide operating range.We conclude that both models can be used for modeling and designing MOSFET integrated circuits for operating frequencies of at least 1.2 THz (of special interest for 6G communications, THz imaging, biomedical diagnostics, etc.).Finally, it is worth to mention that the ADS-HDM can be adapted to several key material technologies in the THz research field such as graphene, AlGaN/GaN-(as shown in [21]) and AlGaAs/GaAs HEMTs, as the governing hydrodynamic transport equations only need to be coupled to an charge control model suitable for the respective material technology.

A Lock-in amplifier correction factor -Toptica TeraScan 1550
In Fig. 6(a), we sketch the measurement setup used to gather the data in this work.The Si CMOS TeraFETs are characterized using the lock-in technique (Ametek DSP 7265) with an electronically modulated InGaAs photomixer emitter (Toptica TeraScan 1550) @ f mod = 7.62 kHz (indicated as blue line).Alignment of the detectors was performed in xyz to guarantee the correct position in the focal point of the final OAP mirror (diameter: 2").From time-domain acquisitions of the detector signal (red line), we determined the lock-in correction factor α LIA via  as a function of the radiation frequency using a single broadband 65-nm Si CMOS TeraFET (log-spiral antenna).In Fig. 6(b), the measured detector signal U DS (red line, left axis), the modulation waveform for the Toptica emitter (blue line, right axis) -which serves as the reference for the lock-in amplifier -, and the fundamental component of the measured detector signal U lock−in,1f mod (green line, left axis) -extracted using Fourier analysis -are presented for three different radiation frequencies (0.2, 0.6 and 1 THz).In addition, the extracted peak-to-peak value (red dashed line), relevant for the detector calibration and sensitivity measurements, as well as the obtained RMS value of the fundamental component (green dashed line), which is measured by the lock-in amplifier after bandpass filtering at f mod = 7.62 kHz, are shown.In Fig. 6(c), we present α LIA,corr as a function of the radiation frequency.A linear fit is used to determine the observed frequency trend, which is relevant for our simulations and correct predictions of ∆U DS .We speculate that the unexpected frequency dependence of α LIA,corr is a consequence of frequency-dependent non-linear mixing properties of the Toptica InGaAs photomixer emitter.They manifest, when the photomixer is illuminated by the fiber-coupled dual-color optical beam at different beat-note frequencies ∆ν (0.1-1.2 THz).This leads to the observed non-ideal (here non-sinusoidal) modulation of the emitted THz output power and to the change in α LIA,corr with the beat-note frequency.Further investigations, also at higher frequencies, are needed to clarify the origin of the observed effect.

B Gaussian beam coupling efficiency
In order to calculate the total power coupled into the effective area (indicated as green square in Fig. 7) of an arbitrary planar antenna we assume a point symmetry of I gauss (x, y) to split the calculation into 8 equal triangles as indicated in

Figure 1 :
Figure1: (a) Schematic of the ADS-HDM simulation tool: The blue box (full simulation domain) represents the THz detector model (is taken in the same manner for the TSMC RF foundry model), here for source radiation coupling conditions.At its core (hydrodynamic domain), the hydrodynamic transport equations are used to calculate the rectified signal of the detector on the basis of plasma-wave-modified resistive mixing in the transistor's channel.The equations are solved in the Keysight ADS simulator itself.The FET parasitics are taken into account in the FET circuit domain.Parasitic components considered are the source and drain overlap C ov and fringe capacitances C f r , the contact resistances R C,S and R C,D , and the gate resistance R G .The antenna parameters including the antenna impedance Z ant are estimated from EM wave simulations (CST Studio Suite).For details of the HDM implementation, we refer to[21].(b) Cross-sectional view of the detector implementation in TSMC 65-nm process.The patch antenna is implemented underneath a passivation layer in the top metal layer MT.The ground plane is placed in lowest to metal layers M1 and M2.

Figure 2 :
Figure 2: (a) Top-view micrograph of the array of four detectors (A1-A4) showing the rectangular patch antennas; the corresponding patch lengths are listed.All patches have the same width of 40 µm and height of 8.41 µm above the ground-plane.(b) Direct measurement of the Gaussian beam intensity profile I gauss (x, y) obtained close to the resonant frequencies of the patch antennas A1-A4.The physical dimensions (L patch • W patch ) of each patch are indicated as blue rectangles, while the simulated effective area A ef f is shown in green.(c) Simulated real (solid line) and imaginary part (dashed line) of the antenna impedance Z ant and (d) effective area A ef f obtained from EM wave simulations (CST Studio Suite).(e) Extracted FWHM in x (red), y (blue) and estimated mean FWHM (green) of the Gaussian beam intensity profile measured in the focal point of a 2" OAP mirror as presented in (b).The black solid line shows the calculated minimum FWHM from Eq. (7) for d L = 30 mm and f L = 50.8mm.(f) Antenna Gaussian beam coupling efficiency η gauss for each detector (A1-A4) extracted from direct measurements of the Gaussian beam intensity profile (open dots).In addition we present η gauss for arbitrary frequencies assuming (i) an ideal Gaussian beam intensity profile (open squares) and (ii) solving Eq. (9) numerically (solid lines).In both cases the FWHM of the Gaussian beam is determined from Eq. (7) (as shown as black line in (e)).

Figure 3 :
Figure 3: Comparison of the measured DC drain-to-source resistance (black line) with simulations obtained from TSMC RF (blue open squares) and ADS-HDM (red open squares) for the different Si CMOS TeraFETs (A1-A4).U DS is set to 10 mV for all measurements and simulations.The dashed vertical line represents the extracted U T for the respective devices using a fit model based on the UCCM (see Table1).For both models we assume a fixed residual resistance of R p,ext = 300 Ω for a short-circuit protection FET, which is connected in series with the drain electrode of all respective transistors.

Figure 4 :
Figure 4: Top panels: Comparison of the measured frequency-dependent optical voltage responsivity (black line) with simulations using TSMC RF (blue open squares) and ADS-HDM (red open squares) for devices A1-A4.Bottom panels:Comparison of the simulated frequency-dependent real (solid lines) and imaginary (dashed lines) antenna impedance Z ant (green) with corresponding simulations of the complex conjugate transistor channel impedance Z t * using TSMC RF (blue) and ADS-HDM (red) for the different Si CMOS TeraFET devices.Estimated impedance mismatch factor η m (dots, right axis) is calculated from Eq.(10) for TSMC RF (blue) and ADS-HDM (red).

12 Figure 5 :
Figure 5: Left panels: (a) Comparison of the measured voltage-dependent optical responsivity R op V (black line) at constant excitation frequency with corresponding simulations obtained from TSMC RF (blue open squares) and ADS-HDM (red open squares) for A1-A4.The vertical dashed line represents the extracted U T from R DS measurements in case of the ADS-HDM (see Table 1).The pure resistive self-mixing detector response predicted by the ADS-HDM is plotted as red line.Right panels: (b) Comparison of the measured optical voltage-dependent NEP (NEP op V ) (black line) at constant excitation frequency with corresponding simulations using the TSMC RF foundry model (blue open squares) and our in-house tool ADS-HDM (red open squares) for the different Si CMOS TeraFET device.The vertical dashed-dotted line represents minimum NEP opV of the respective devices.TSMC RF (blue filled squares) and ADS-HDM (red filled squares) simulations of the voltage-dependent electrical NEP (NEP el V ) with regards to the maximum available power delivered from the antenna P el in,ant .In addition, the simulated electrical NEP calculated with regards to the available power at the source transistor terminal P el in,t is depicted for TSMC RF (blue dashed line) and ADS-HDM (red dashed line).
Figure 5: Left panels: (a) Comparison of the measured voltage-dependent optical responsivity R op V (black line) at constant excitation frequency with corresponding simulations obtained from TSMC RF (blue open squares) and ADS-HDM (red open squares) for A1-A4.The vertical dashed line represents the extracted U T from R DS measurements in case of the ADS-HDM (see Table 1).The pure resistive self-mixing detector response predicted by the ADS-HDM is plotted as red line.Right panels: (b) Comparison of the measured optical voltage-dependent NEP (NEP op V ) (black line) at constant excitation frequency with corresponding simulations using the TSMC RF foundry model (blue open squares) and our in-house tool ADS-HDM (red open squares) for the different Si CMOS TeraFET device.The vertical dashed-dotted line represents minimum NEP opV of the respective devices.TSMC RF (blue filled squares) and ADS-HDM (red filled squares) simulations of the voltage-dependent electrical NEP (NEP el V ) with regards to the maximum available power delivered from the antenna P el in,ant .In addition, the simulated electrical NEP calculated with regards to the available power at the source transistor terminal P el in,t is depicted for TSMC RF (blue dashed line) and ADS-HDM (red dashed line).

Figure 6 :
Figure 6: (a) Schematic drawing of the measurement setup for optical characterization of the Si CMOS TeraFETs.(b) Fundamental component (green) of the detected waveform (red) extracted from time-domain acquisitions (using an oscilloscope) for a single broadband 65-nm Si CMOS TeraFET at 0.2, 0.6 and 1 THz when illuminating the detector with an electronically chopped (reference shown in blue) InGaAs photomixer (Toptica TeraScan 1550).(c) Extracted lock-in correction factor α LIA,corr as a function of frequency.Calculated values (dashed / dotted lines) for the lock-in correction factor are shown for different types of references/emitter modulations.A reference measurement of α LIA,corr at 1 THz measured with the similar Toptica TeraScan 1550 is shown as blue star[39].

( 13 )LFigure 7 :
Figure 7: Sketch for visualizing mathematical steps for the calculation of the Gaussian beam coupling efficiency.The dimensions of the effective area A ef f = L 2 are indicated as a green square.

Table 1 :
Extracted simulation parameters for the ADS-HDM of the different Si CMOS TeraFETs obtained through a fitting model based on the UCCM.W is the width of the channel.R C is the contact resistance for Source and Drain contact respectively.γ is the non-ideality factor.

Table 2 :
Measured minimum NEP op V in pW/√Hz at respective gate voltages U GS for the detectors A1-A4 as well as simulated values obtained from TSMC RF and ADS-HDM for NEP op V (op.) and NEP el V (el.) with regards to P el in,ant .op = 1, η gauss = 1, η scat = 1) and 10 − 12 pW/