We propose a graphene-based terahertz detection scheme capable of measuring not only pulse energy but also electric field shape. The scheme leverages strong nonlinear velocity saturation characteristics of graphene in combination with envelope-carrier phase offset imposed by propagation of the pulses through the dispersive medium to produce shape-dependent electric charge. These charges can then be easily detected by conventional electronics, and as numerical calculations show, the original pulse shape can be recovered with the help of deep neural networks.

Terahertz (THz) technology is gaining ever wider acceptance in applications such as material characterization,^{1} study of carrier dynamics in nanostructures,^{2} spectroscopy,^{3} and imaging.^{4,5} THz spectroscopy is especially interesting for biological applications, due to its noninvasive and nondestructive nature,^{6} as well as for applications such as detection of explosives.^{7} Unfortunately, further advances in THz spectroscopy are impeded by a lack of efficient detectors^{8} and powerful sources. Since most of the THz sources are pulsed, THz measurements are performed in time domain with subsequent spectral analysis performed digitally, but typical THz detectors rely on slow pyroelectric or bolometric response^{9} and often require cooling to cryogenic temperatures. In the absence of fast THz detectors for time domain measurements, THz signals are usually upconverted via interactions in nonlinear crystals and subsequently detected by optical detectors.^{10} Alternatively, one can use photoconductive switches,^{11} detectors based on plasmonics,^{12} or the Dyakonov-Shur effect.^{13} Since neither of these detectors is capable of measuring short (subpicosecond) time intervals, optical delay lines are employed to measure the time indirectly with high resolution which greatly increases the complexity of THz measurements. THz photodetectors based on graphene and other two dimensional materials have demonstrated a wide bandwidth, excellent sensitivity, and fast response^{14–20} as discussed in detail in a superb review by Koppens *et al.*^{21} However, they are still quadratic detectors and can only provide information about the field energy.

In this work, we propose a fundamentally different THz detection scheme. In this scheme, a drift current is induced by the THz electric field in a graphene monolayer (or any other semiconductor with high mobility and substantial velocity saturation) at room temperature. This current is then integrated by a relatively slow electronic circuit, and the resulting charge is measured. This is analogous to measuring ultrafast current pulses produced by a different phenomenon—coherent photoinjection in Refs. 22 and 23. The temporal field profile is recovered by reconstruction from a series of distinct measurements where the THz pulse is propagated through a highly dispersive medium. Such propagation conditions cause changes to the pulse profile, and as a result, it leads to changes in the detected THz-induced charge.

The proposed detection apparatus is shown in Fig. 1 and consists of a wedge with variable thickness *L* moving in and out of the THz beam path, which introduces variable dispersion, followed by a graphene detector. Of course, one can also employ a graphene detector array if a single-shot measurement is desired. The detector itself, shown in the inset, is a simple graphene square of size *W *×* W* commensurate with the diffraction spot size of the focusing system (i.e., a few hundred micrometers) placed on the insulating THz-transparent substrate. The detector is connected to an amplifying circuit with input capacitance *C _{in}* via two strips making Ohmic contact with graphene. Graphene itself is conductive and doped to two-dimensional carrier density n (possibly electrostatically via the backgating through the substrate).

The impinging THz field *E*(*t*) generates a current of free carriers in the graphene monolayer

where $\mu (E)$ is the mobility and *q* is the electron charge which is the dominant source of photocurrent in the regime $\u210f\omega <kBT$ (where electric field shapes can be routinely measured as done in electronics) with $\u210f$ being the Planck constant, *ω* the electric field frequency, *k _{B}* the Boltzmann constant, and

*T*is the temperature. The graphene strip acts as an AC source; however, this current is difficult if not impossible to detect as the signal gets averaged by an amplifier having feedback capacitance

*C*and input resistance

*R*. The generated high frequency current gets convolved with the detector response

where the time constant $\tau =RC$ is orders of magnitude longer than the THz period *T*. Consequently, for the few cycle THz pulse of length *T _{p}*, what is being measured is essentially a charge

If the velocity-field characteristic is linear (i.e., mobility is field independent), the integral vanishes since in the optical (THz) domain, the average electric field is always 0, $\u222bE(t)dt=0$. Even if the average THz polarization at the THz source is not zero, the electric field at low frequencies diffracts away before it can reach the detector. If, on the other hand, the velocity field dependence is nonlinear, i.e., the mobility is field dependent, the integral does not vanish because in general $\u222bE\alpha (t)dt\u22600$. For example, a detector with *α* = 2 (a typical square law detector such as diode) provides information about the pulse energy *U _{p}* as

A square law detector, however, can only provide information about the pulse energy. A higher even-order $\alpha =2m$ response would add information about the temporal extent of the pulse but not about its shape. However, for odd $\alpha =2m+1$, far more information about the pulse temporal shape can be extracted as shown in Refs. 24 and 25 for *α* = 3 and in Ref. 26 for higher odd *α*'s. For a signal described by the product of a THz carrier and envelope $E(t)=e(t)\u2009cos(\omega 0t+\delta \varphi )$, the value of integral $\u222bE2m+1dt$ depends on the relative carrier-envelope phase offset $\delta \varphi $. Assuming the pulse envelope *e*(*t*) to be symmetric, the integral is equal to zero whenever $\delta \varphi =(k+1/2)\pi $ and reaches its maxima and minima whenever $\delta \varphi =k\pi $.

The change in phase offset can be implemented by propagating the pulse through the dispersive medium of length *L*. In the dispersive medium, the phase of the carrier propagates with the phase velocity $vp=\omega /k$ while the peak of the envelope moves with the group velocity $vg=d\omega /dk$ so that $\delta \varphi =\omega L(vp\u22121\u2212vg\u22121)$ which leads to an oscillating detected voltage as a function of *L*.

A medium with strong odd-order nonlinear response can yield distinct values for different field profiles which should be enough to detect high intensity THz field shapes. To that end, the obvious place to look is a material with strong velocity saturation, which usually can be traced to band nonparabolicity and therefore to narrow bandgap semiconductors. Naturally, the ultimate material with the strongest velocity saturation is graphene. Measurements of drift velocity for graphene have been described well experimentally using the following expression:^{27}

where *γ* is a fitting parameter ranging between 0.7 and 1.2 and *v _{sat}* is the saturation velocity. Only the field strength interacting with the graphene free electrons is considered here, which is a fraction of the incident field. As shown in the supplementary material, for our example, the incident electric field is 3.5 times higher than the effective field. Although absorbed power is quite significant (1 W), it is not sufficient to melt either the substrate or the graphene layer (see the supplementary material). For maximal nonlinearity, we assume $\gamma =1.2$ which yields the current-field characteristics shown in Fig. 1(b) for the saturation velocity of $vsat=7.8\xd7107\u2009cm/s$, the saturation field strength of $Es=1\u2009kV/cm,$ the device width of $W=100\u2009\u2009\mu m$, and a carrier density of $nd=0.8\xd71012\u2009cm\u22122$ which results in a saturation current of $Is=nd\xd7vsat\xd7W\xd7q\u2248100\u2009\u2009mA$. The THz source must therefore have field strengths larger than $3.5\u2009\u2009kV/cm$ which can be attained experimentally,

^{28}and theoretically, even larger fields should be possible.

^{29}Polyethylene naphthalate (PEN), whose refractive index $n(\omega )$ is shown Fig. 1(c), can be used as the dispersive material since it is transparent in the range of interest.

^{30}Suppose we would like to calculate the detected signal for the following pulse $E(t)=Ae\u2212(t/tw)2\u2009sin(\omega 0t)$, where

*A*is the field amplitude,

*t*is the temporal width, and

_{w}*ω*

_{0}is the carrier frequency. The Fourier spectrum of the pulse is $E(\omega )$, and after propagating distance

*z*through the PEN prism, an additional phase $\varphi (\omega )=\omega cn\u0303(\omega )z$ is acquired. Hence, the time domain shape of pulse can be found as

The evolution of pulse shape is shown in Fig. 2(a) for three propagation distances corresponding to the different carrier-envelope phase offsets. As one can see, the initially perfectly odd shape on the left becomes even when phase offset is equal to 90° and then once again becomes odd at twice the distance (180° phase offset).

The corresponding current pulses are shown below the pulse shapes in Fig. 2(b), and the velocity saturation is clearly seen as the current amplitudes get suppressed. At the same time, the symmetry of the pulses is preserved, and hence as a result of integration shown on the line below [Figs. 2(c)–2(d)], the accumulated charge *Q* is zero for the odd shape and nonzero for the even shape (although not exactly zero for the third column due to higher order dispersion). Therefore, as one plots the accumulated charge (or voltage on the capacitor $V=Q/C$) as a function of the propagation length *L* [Fig. 2(e)], an oscillatory response is obtained. The oscillation period is related to the phase and group velocity difference at a given carrier frequency *ω*_{0}. The amplitude of the oscillations decays due to dispersion-induced pulse broadening which depends on the group velocity dispersion (GVD) and the pulse width. The prism length in this proof of concept numerical demonstration is large since the refractive index profile is very flat. A flat profile is picked for lower loss and faster numerical calculations of field shape. In practice, materials (or diffraction gratings or even metasurfaces) with larger dispersion can be used to keep the disperser dimensions small.

Signal levels from the detector is of the order of microvolts $V\u223c9.3\xd710\u221218C10\u221212\u2009F$ as seen in Fig. 2 for a detection capacitance of 1 pF, an amplitude of $1.5\u2009kV/cm$, a center frequency of 2 THz, and a temporal width of 0.7 ps. Calculations in the supplementary material show that voltage persists for roughly 100 ps, and hence, electronics operating at 10 GHz are needed. The detected signal is not always sinusoidal due to higher order dispersion. It is easy to see that higher Taylor components in the dispersion relation provide additional information about the pulse shape. For example, group velocity dispersion (GVD) leads to pulse broadening, and if the pulse is sufficiently broadened, the carrier drift velocity will approach the linear regime at some point yielding zero integrated charge. Therefore, the integrated charge will decay to zero with a characteristic decay rate determined by GVD which is known, and the pulse width can be recovered.

The complex relationships between the input pulse shape, material dispersion, and integrated charge output cannot be expressed analytically. Therefore, one must resort to numerical methods to solve the inverse problem. One solution would be to define some metric such as mean squared error between a measurement and a model prediction, which can then be minimized using model parameters. However, running the optimization algorithm would take prohibitively long, making the detector infeasible. Instead, we propose using deep neural networks (DNNs), which can carry out inference in real time, since the functional dependence between the integrated charge and the pulse shape can be learned with sufficient training data. As a proof of concept, we simulated the field detection scheme by taking the first two Taylor components of the dispersion relationship (i.e., phase and group velocity) for PEN, which is a reasonable approximation in the frequency range of interest (2–3 THz), and trained a DNN on the measured voltages. DNNs require a lot of data for training, which can be generated very quickly numerically if only the phase and group velocities are considered.

The neural network shown in Fig. 3 is trained using PyTorch, a popular open source machine learning library, on a commercially available GPU (graphics processing unit) with mean squared error (MSE) as the loss function. Pytorch implementation of ASGD (averaged stochastic gradient descent)^{31} with the default values (a learning rate of $10\u22122,\u2009\lambda =10\u22124,\u2009\alpha =0.75$ and weight decay parameter of 0) is used for optimization. The training set consisted of 100 000 samples drawn randomly from a uniform distribution of center frequencies in the range of 2–3 THz. The electric field amplitude is varied between 3.5 and $4.5\u2009\u2009kV/cm$, and the temporal width is changed so as to keep the total interacting power constant ($tw[ps]\xd7A2[kV2/cm2]=16$). The input layer of the DNN takes in voltage measurements such as the one showed in Fig. 2(e) for $L=10\u2009cm$ discretized in space (100 data points). We found that the problem of the inferring pulse width and the center frequency is ill-posed when only a single detector is used, and a significant improvement in the reconstruction is observed when another set of measurements from a detector with a different saturation curve becomes available. This can be implemented, for example, by changing the carrier concentration in the graphene sheet. Increasing carrier density will reduce the saturation field strength (by half in our example) and increase the saturation current (by a factor of two in this example). This can be experimentally achieved by back gating the same graphene monolayer. Measurements from both sets are then concatenated and input to the DNN.

The DNN prediction and the ground truth (original pulse) have excellent agreement as shown in Fig. 4(a) for a center frequency of 2 THz, a pulse width of 1 ps, and an electric field amplitude of $4\u2009\u2009kV/cm$. In order to study the DNN performance more rigorously, 1000 test samples were drawn randomly from the same distribution, and reconstruction errors are calculated defined as $f0,gnd\u2212f0,pre$ and $tw,gnd\u2212tw,pre$ for the center frequency and the pulse width, respectively, where the subscripts “gnd” and “pre” correspond to ground truth and predicted values, respectively. Histograms of errors are shown in Figs. 4(b) and 4(c) with mean absolute errors of 0.067 THz and 0.041 ps, respectively. Although we have only used a Gaussian function, it is possible to train the network with a combination of any suitable basis functions such as Lorentzian and Voigt profiles. The combination that best describes a particular THz source does not need to be known *a priori* since layers similar to the inception blocks in GoogLeNet can be used,^{32} where all suitable basis functions are considered at once. During training, the network would then automatically pick the set that best describes the source. The network could also be trained with reference-sample arms to perform spectroscopic measurements.

In conclusion, we have shown in this proof of concept numerical study that DNNs can be used to infer the ultrafast field shapes from a series of distinct measurements, relying on nonlinear response of graphene and controlled distortion of the field using dispersion. It should be noted that the DNNs do not require any prior information about the dispersion curve, the thickness of the prism, or any other physical parameter. Therefore, an experimental implementation does not necessarily require any knowledge of the physical parameters but simply some ground truth input-output pairs for training. However, in many cases, generation of training data experimentally might be challenging, in which case the training data can be generated numerically given physical parameters and models as done here. Once trained, the network can perform the inference in real time.

See the supplementary material for detailed calculations of reflection, transmission, and absorption coefficients; electron and lattice temperature change due to absorption; and Python implementation of the deep neural network.

This work was funded by NSF Award No. 1741694.