Terahertz Nonlinear Optics of Graphene: From Saturable Absorption to High‐Harmonics Generation

Graphene has long been predicted to show exceptional nonlinear optical properties, especially in the technologically important terahertz (THz) frequency range. Recent experiments have shown that this atomically thin material indeed exhibits possibly the largest nonlinear coefficients of any material known to date, paving the way for practical graphene‐based applications in ultrafast (opto‐)electronics operating at THz rates. Here the advances in the booming field of nonlinear THz optics of graphene are reported, and the state‐of‐the‐art understanding of the nature of the nonlinear interaction of graphene with the THz fields based on the thermodynamic model of electron transport in graphene is described. A comparison between different mechanisms of nonlinear interaction of graphene with light fields in THz, infrared, and visible frequency ranges is also provided. Finally, the perspectives for the expected technological applications of graphene based on its extraordinary THz nonlinear properties are summarized. This report covers the evolution of the field of THz nonlinear optics of graphene from the very pioneering to the state‐of‐the‐art works. It also serves as a concise overview of the current understanding of THz nonlinear optics of graphene and as a compact reference for researchers entering the field, as well as for the technology developers.

that is, the multi plication of the frequency of incident THz photons in graphene, remained elusive. This situation has changed recently when very efficient generation of THz harmonics up to the seventh order in single-layer-doped graphene was demonstrated using multicycle quasi-monochromatic THz driving fields, revealing gigantic nonlinear optical coefficients of graphene in the THz range. [63] The physical mechanism of the demonstrated highly efficient THz high-harmonics generation (THz HHG) and THz saturable absorption in doped graphene was attributed to the collective thermodynamic response of the background Dirac electrons of the doped graphene to the THz driving field, which is currently believed to be the dominating effect in interaction of free electrons in graphene with THz fields. The underlying physical picture of this interaction, and its effect on the THz conductivity and absorbance in graphene will be provided within this report. Finally, we emphasize the importance of exploring the THz nonlinear properties of graphene for emerging technologies. Graphene-based transistors, [64][65][66][67] photodetectors, [68,69] saturable absorbers for laser modelocking, [30,70,71] and other high-speed (opto-)electronic devices, all usually operate in the nonlinear regime, and are entering the THz frequency/rate operation mode.

How Is This Report Organized?
Our report is organized in five sections. After the introduction in section 1, section 2 provides a brief overview of the important fundamentals of graphene, such as its band structure and key electronic properties. Section 3 deals with the linear interaction of graphene with light fields in general, accounting for the interband and intraband electron dynamics. Section 4 presents (i) the experimental observations of nonlinear THz optics of graphene and other graphene-based systems, and a brief overview of the experimental techniques and facilities for the generation of intense THz fields, (ii) the thermodynamic model that describes the nonlinear interaction of graphene with intense THz fields, and (iii) the experimental generation of THz high harmonics, which directly follows from the thermodynamic picture of THz electron transport in graphene. Where relevant, the results on graphene under optical excitation will also be briefly discussed. Finally, section 5 provides conclusions and future perspectives on possible applications of graphene as a highly nonlinear THz material.

Band Structure of Graphene
Graphene is a monolayer of carbon atoms bound covalently in a hexagonal honeycomb lattice, as shown in Figure 1a. Each carbon atom in the graphene plane is bound to three nearest neighbor atoms via strong equivalent σ-bonds established by sp 2 hybridization of three of the valence electron-orbitals of each atom, namely, 2s, 2p x , and 2p y , [26][27][28]72] as shown in Figure 1b. These in-plane σ-bonds are responsible for the mechanical strength of graphene. The fourth electron-orbital 2p z is oriented perpendicular to the graphene plane. The overlap between the out-of-plane 2p z orbitals of the carbon atoms results in weaker π-bonds. These delocalized π-electrons are responsible for the electrical conduction in graphene.
The hexagonal lattice of graphene is considered as a bipartite lattice composed of two triangular sublattices, marked A and B in Figure 1c. The sublattices are the Bravais lattices of graphene with lattice vectors = 0.5 ( 3, 1) 1 a a and = − 0.5 ( 3, 1) 2 a a and a lattice constant a = |a 1 | = |a 2 | = 0.246 nm. [27] We note here that the lattice constant a is related to the carbon-carbon distance a cc = 0.142 nm by = 3 cc a a . In the reciprocal space, the first Brillouin zone is consequently a hexagon featuring six corner points of two inequivalent groups of K-points called the Dirac points indicative of the two sublattices A and B, as shown in Figure 1d. The corresponding band structure is shown in Figure 1e, in which the valence band touches the conduction band at these six Dirac points, providing a gapless band structure of graphene. This band structure of graphene was first calculated by Wallace in 1947 using the tight-binding approximation, [73] providing cosine-like energy bands in energy-momentum space [27] k k k a k a k a x y where γ 0 ≈ 2.8 eV is the hopping energy (transfer integral) between the nearest-neighbor atoms. Near the Dirac points where |k|a ≪ 1, the Taylor series expansion of the energy spectrum of Equation (1) reduces to the linear dispersion relation [27]  ε γ where  γ = ≈ 3 /2 10 F 0 6 v a m s −1 is the energy-and momentumindependent band velocity (Fermi velocity) [73] and the ± sign is indicative of the positive kinetic energy of electrons in the conduction band (the π* electrons) and the negative (binding) energy of the bound electrons (the π electrons) in the valence band. This unique linear dispersion ε = v F ℏk with constant band velocity in the vicinity of the Dirac points in graphene is characteristic of ultrarelativistic massless quasiparticles: Dirac fermions with zero rest mass. The linear band structure approximation of Equation (2) is valid for the energies up to about ± 1.5 eV with respect to the Dirac point, thus covering a very broad range of energies (frequencies) of the electromagnetic spectrum from THz to UV.
The density of states in graphene D(ε) within the linear part of the band structure is a linear function of energy, with D = 0 at the Dirac point ε = 0  ε π ε ( ) = The Fermi energy E F depends on the free carrier density N c in graphene as Very importantly for the electronic transport in graphene, it can be shown that for the "massless" electrons with linear energy-momentum dispersion such as Equation (2) it is fundamentally not possible to electrostatically bind them to a particle of opposite charge. This is in stark contrast to the "massive" charge carriers with quasi-parabolic energy-momentum dispersion, such as electrons and holes in regular metals or semiconductors, which can be trapped by charged impurities, or form excitons. This absence of electrostatic (Coulomb) binding for its electrons also makes graphene stand out among other 2D materials with "massive" quasi-parabolic energy momentum dispersion for the charge carriers, such as transition metal dichalcogenides, which feature extremely strong Figure 1. Graphene lattice and band structure: a) Hexagonal honeycomb lattice of graphene, b) binding between the carbon atoms in the graphene lattice, c) the triangular sublattices A and B, d) the first hexagonal Brillouin zone exhibiting two inequivalent groups of Dirac points, and e) the energymomentum spectrum of the graphene band structure described by Equation (1). Note the rotation of k-axes with respect to (d). (Reproduced, except for (b), with permission. [26] Copyright 2009, American Physical Society). electronic Coulomb binding due to the very weak dielectric (lattice) screening, leading in some cases to exciton binding energies as high as ≈500 meV. [74] In graphene, due to the fundamental absence of Coulomb binding, electron transport is not affected by the trapping onto impurities and other unscreened Coulomb potentials, and is only limited by the momentum scattering. This unique property of graphene electrons, combined with the large optical phonon energy of graphene of the order of ≈200 meV, and the nonpolar nature of its lattice, leads to a very large carrier mobility that can exceed 100 000 cm 2 V −1 s −1 observable even at room temperature. [66,72,75] Further, graphene exhibits a strong ambipolar field effect when subjected to a gating voltage, allowing for a precise control of the type and density of the free charge carriers (electrons or holes) in the respective bands. [1] Combined with naturally high electron mobility, this made graphene a promising candidate for the development of ultrafast transistors operating at the rates of 100s of GHz. [64,65]

Light-Matter Interaction in Graphene in the Linear Regime
Depending on the photon energy ℏω and the Fermi energy of graphene E F , the excitation of graphene with light will lead to either interband or intraband transitions. The pioneering observation and distinction between both types of transitions in graphene was made by Winnerl et al. [9] Optical excitation of graphene with ℏω > 2E F results in an interband transition, similar to photoexcitation of semiconductors, and leads to the creation of an electron-hole pair, as shown schematically in Figure 2a.
The interband transitions in graphene exhibit a universal, that is, frequency-independent, absorption A = πα = 2.3% of the incident light, [5] where α = (1/4πε 0 ) e 2 /ℏc = 1/137 is the finestructure constant, corresponding to the universal optical conductivity σ 0 = e 2 /4ℏ  1 G. Such a universal light absorption is a direct consequence of the unique electronic properties of graphene: the linear energy-momentum dispersion and the gapless band structure. We note that for photon energies below ≈0.5 eV and above ≈2 eV, a deviation from this universal absorption occurs. [6,12,76] The frequency dependence of the corresponding interband conductivity was shown to follow the relation [6,[8][9][10]12] where k B is the Boltzmann constant and T e is the electron temperature. We note that for large photon energies, that is, ℏω ≫ 2E F and ℏω ≫ k B T e , Equation (5) reduces to the universal conductivity σ 0 = e 2 /4ℏ. We further note that the optical absorption (and reflection) of graphene depends strongly on its dielectric environment. [77,78] For photons of energy ℏω < 2E F , as is typically the case at THz frequencies, the interband transitions are Pauliblocked, [9,12,[79][80][81] whereas the intraband transitions are possible, as depicted schematically in Figure 2b. Such transitions correspond to the interaction of the incident light with the electrons around the Fermi surface of graphene, that is, to the free carrier absorption. The frequency-dependent optical conductivity corresponding to such intraband transitions, in case of a thermally equilibrated electron population, is described by the following solution of the Boltzmann transport equation [28,59] where τ(ε) is the energy-dependent electron momentum scattering time, and f FD (ε, μ, T e ) is Fermi-Dirac distribution function depending on the electron energy ε, chemical potential μ, and electron temperature T e . For the case of T e ≪ E F /k B , which corresponds to a sharp Fermi edge, only the electron scattering time at the Fermi energy, τ = τ(E F ), matters, and Equation (6) reduces to the following Drude-type expression [8,9,82] σ The THz excitation of graphene is dominated by the intraband conductivity (Equations (6) and (7)), especially at room temperature. Indeed, the photon energy in the few meV range (1 THz = 4.1 meV) is considerably smaller than the roomtemperature thermal energy k B T e = 26 meV. As a result, in intrinsic, undoped graphene with E F ≈ 0, at room temperature, the thermal population of the conduction and valence bands becomes almost equal: (f FD, conduction − f FD,valence ) | 1THz, 300K = 0.04 for the interband transition at the frequency of 1 THz = 4.1 meV. This leads to effective Pauli-blocking of the interband THz transitions in graphene at room temperature, making them only possible under cryogenic conditions (see, e.g., ref. [62]). However, in practical THz experiments on large area, ≈mm sized, single-layer graphene samples, even at lower temperatures the THz interband transitions are not easy to reveal because of the presence of so called "charge puddles," local areas of n-or p-doping with |E F |> 0, again Pauli-blocking the interband THz transitions. [83] For the typical doping levels for substrate-supported graphene, resulting in a Fermi energy on the order of 100 meV (see, e.g., refs. [59,63,72,[84][85][86][87][88][89][90]), the interband THz transitions are Pauli-blocked even at T = 0, leaving only the intraband transitions described by the free carrier conductivity of Equation (6) possible.
We note that most single-layer graphene samples are in fact unintentionally doped due to impurities induced in the preparation process as well as the effects of the environment in which graphene is used. [72,[85][86][87][88][89][90] For example, in the course of preparation of chemical vapor deposition (CVD) grown graphene, which is one of the most popular graphene forms broadly available nowadays, hole-doping is induced during the transfer of poly(methyl methacrylate) (PMMA)-covered graphene from copper foils to dielectric substrates by residues of PMMA. [90] The dielectric substrate itself, and water and oxygen molecules trapped between the graphene and the substrate further contribute to the hole-doping of CVD graphene. [90,91] On the other hand, in epitaxially grown graphene on SiC electrondoping is induced by the SiC substrate (substrate-induced cha rging). [9,27,86,92,93] For multilayer epitaxial graphene, the layers closer to the middle of the sample can be intrinsic with E F ≈ 0, while outer layers are usually doped by the environment and the substrate. [9,62,93] Given the symmetry of Dirac band structure of graphene (see Figure 1), there is no fundamental difference between the electron and hole doping: Both charge species can be considered equivalent and conduct electric currents in the same fashion. In the literature, for simplicity, the Fermi level is often indicated within the conduction band even if graphene is in fact hole-doped (such as typical CVD-grown graphene). Correspondingly, irrespective of the actual type of majority carriers in graphene, they are often referred to as electrons. In this paper, we follow this convention, unless specifically indicated.
Figure 2c-f shows the THz field transmission through a monolayer CVD-grown p-doped graphene sample, and the corresponding THz power transmission and conductivity spectra of graphene. With THz time-domain spectroscopy (THz-TDS), [94,95] which is based on direct time-domain sampling of the propagating THz electric field, one can extract the graphene conductivity from the experimentally measured THz field transmission of the graphene film through the following relation in the frequency domain (known as Tinkham equation) [  Interaction of graphene with a light field: Schematic of the band structure with a) interband transition due to absorption of a photon of energy larger than twice the Fermi energy, ℏω > 2E F , leading to electron-hole pair generation, followed by the energy and momentum exchange with the Fermi sea electrons via electron-electron scattering, and b) intraband dynamics when ℏω < 2E F , which exhibits Drude-like conductance. c) Demonstration of the transmission of a THz field through a graphene sample with absorbed portion of the field ΔE THz related to the graphene intraband conductivity (Reproduced under the terms of the Creative Commons CC-BY 4.0 License. [59] Copyright 2015, Springer Nature). d,e) Actual measurement of THz fields and the corresponding spectra, respectively, transmitted through a graphene sample on SiO 2 substrate and through a bare SiO 2 substrate as a reference, and f) the corresponding complex THz intraband conductivity normalized to the universal conductivity σ 0 = e 2 /4ℏ, extracted from the complex spectra (the Fourier transforms) of the measured fields in (d) using Equation (8) and described by the Drude-like intraband conductivity of Equation (7) using E F = 170 meV and τ = 47 fs as fitting parameters (Reproduced with permission. [63] Copyright 2018, Springer Nature).
measured time-domain fields E s (t) and E ref (t), respectively, as shown in Figure 2d. Here, E s (t) is the THz field transmitted through graphene on a substrate, and E ref (t) is that transmitted through a bare substrate of identical thickness used as a reference. Figure 2f shows the complex conductivity, real and imaginary, obtained from the experimental THz fields shown in Figure 2d at room temperature, and fitted by the Drude-like conductivity of Equation (7). The only adjustable parameters in this fit are the Fermi energy E F and the momentum scattering time τ. Thus, THz-TDS allows one to unambiguously determine the key parameters of electron transport in graphene in the linear regime-the number density of conducting charge N c , linked to E F by Equation (4), and the average momentum charge carrier scattering time τ. Note that the sheet conductivity (i.e., conductivity times thickness) in units of S = Ω −1 is used in Equation (8).

General Considerations on THz and Optical/IR Nonlinearities in Graphene
Since the advent of graphene, several theories were put forward that predict a highly nonlinear response of this material, in particular in the THz regime. This field of research was initiated by Mikhailov, [16,17,19] who initially used semiclassical calculations and ballistic transport. Later theories were extended to include quantum mechanical considerations, [18,21,22] and more recent theories of graphene nonlinearities furthermore take into account realistic conditions regarding momentum scattering time and ambient temperature. Interestingly, the effect of the elevated electron temperature on graphene nonlinearities was not considered of significant importance in these earlier works. In 2015, Mics et al. [59] showed that the strong THz nonlinearity in graphene is the result of THz-induced carrier heating, which was followed by several other studies, culminating in the observation of extremely efficient THz high-harmonics generation in graphene in 2018 by Hafez et al. [63] Regarding the visible and near-infrared (NIR) regime, the role of the increased carrier temperature on graphene nonlinearity was first addressed in 2018 by Soavi et al. [46] and Baudisch et al. [97] These studies show that properly taking into account the electron temperature leads to a significant adjustment in the theoretically predicted nonlinearities, yielding excellent agreement between theory and experiment. The nonlinearity in graphene is naturally related to the transport of its electrons in the applied electric field of an incident light wave. In order to understand the fundamental difference between graphene nonlinearities in the THz regime and in the optical or NIR regime, we first summarize the relevant timescales of incident light wave oscillations and of electron dynamics in graphene (see Table 1 and Figure 3).
Adv. Optical Mater. 2020, 8,1900771 [97] The dominant dynamics are based on the coherent photon-electron interactions. Since the induced electron current does not exactly follow the incident electric field, redshifted harmonics are generated. b) Illustration of the nonlinear graphene heat response (red) from incident THz light at 1 THz (blue), after Hafez et al. [63] The dominant dynamics are based on the heating-cooling dynamics of the graphene-electron system, which modulate the graphene intraband conductivity and consequently the induced electron current, leading to re-emission at THz higher harmonics. The lengths of the black lines under the figure panels indicate the fundamental timescales of the relevant processes, scaling with the corresponding figure above them.
We remind that since Coulomb trapping for massless Dirac fermions is fundamentally forbidden, electronic transport in graphene is dominated by momentum scattering events. Graphene electrons undergo momentum scattering by longrange and short-range interactions with charged impurities, disorder, and phonons. [28] In substrate-supported graphene, the former typically dominates, and the momentum scattering time is typically on the order of 10s fs, corresponding to mobilities in the range of 1000-10 000 cm 2 V −1 s −1 . If the energy is supplied to the electron system, this energy is ultimately redistributed among all the carriers by electron-electron scattering, leading to an increase of (collective) electron temperature. Whereas this process depends on several parameters, it typically takes some 10s fs. [98][99][100] Highly excited electrons with the energy exceeding that of the optical phonon energy of ≈0.2 eV can also relax by interaction with strongly coupled optical phonons (SCOPs), which is also believed to occur on the timescale of some 10s fs. [100] For a moderately excited electron system (electron temperature of ≈1000 K), carrier relaxation (i.e., cooling) occurs via coupling to (i) graphene acoustic phonons through disorder-assisted scattering, [101] (ii) graphene-optical phonon interactions, [102] and (iii) substrate phonons through near-field interactions. [103,104] At room temperature, these cooling processes typically take place on a timescale of picoseconds, although longer timescales are predicted to occur in the absence of disorder and a substrate. [105] The nature of the nonlinear response in graphene can be divided into ballistic and diffusive, depending on the timescale of electron transport lifetime (i.e., average momentum scattering time) and the duration of the optical cycles of the driving light wave. Very strong nonlinearity in the electronic motion in graphene in the ballistic regime can be achieved making use of the coherent Bloch oscillations. [106,107] They arise from the acceleration of an electron from the Brillouin zone center to the Brillouin zone edge, with π ∆ = k a , in the electric field of the light wave, in a time faster than the average electron momentum scattering time τ = 10-100 fs, thus making electrons oscillate coherently between Δk and −Δk points in the Brillouin zone (here a = 0.246 nm is the lattice constant of graphene). In the optical or IR ranges, the duration of the optical cycle is on the order of a few femtoseconds, which is faster than the average electron momentum relaxation time in graphene (see Table 1), thus ensuring the ballistic motion of an electron on the timescale of the light wave oscillation. From the acceleration theorem that a minimum electric field on the order of 1-10 MV cm −1 is required to reach the regime of coherent Bloch oscillations in graphene with a typical scattering time of τ = 10-100 fs. Here, ℏk is crystal momentum, t is time, e is elementary charge, and E is electric field strength of the light wave. Such rather high peak fields exceeding 1-10 MV cm −1 are routinely attainable for optical and IR signals generated with amplified laser systems. Indeed, high-harmonics generation in graphene using coherent Bloch oscillations mechanism was recently demonstrated with IR fields of 10s MV cm −1 strength by Yoshikawa et al. [47] and Higuchi et al. [108] We note that coherent Bloch oscillations can be induced at lower field strength in systems with reduced Brillouin zone, such as semiconductor superlattices. [109,110] In contrast to IR or optical light waves, THz waves oscillate on a (sub-)picosecond timescale, considerably longer than the typical electron momentum scattering time in graphene which is on the order of 10-100 fs. Therefore, in the THz range the carrier transport typically occurs in the diffusive regime, where electrons undergo several momentum scattering events within the duration of the THz wave oscillation period (see Table 1), and the microscopic motion of electrons is generally incoherent with the driving electric field. As we will explain in detail, the nonlinearity of graphene related to such diffusive collective motion of free electrons in the applied THz fields is in fact extremely strong, leading to substantial THz saturable absorption [9,55,57,59] and extremely efficient THz high-harmonics generation. [63] Remarkably, this nonlinearity mechanism only requires THz fields of the order of 10s of kV cm −1 , which is ≈3 orders of magnitude smaller (i.e., requiring about one million times less peak power) than that required for the optical or IR nonlinearity in graphene via coherent Bloch oscillations. [47,108] In general, THz nonlinearities in graphene can originate from both intraband and interband THz transitions. Simple considerations, however, point out that the THz intraband nonlinear response is expected to be much stronger than the interband one. First, the number of states available for the interband THz transitions is fairly small. Integration of the density of states in graphene (Equation (3)) within 1 THz = 4.1 meV band around the Dirac point yields only about N ≈ 10 8 cm −2 electrons available for interband THz transitions. On the other hand, a typical doping of graphene readily yields a Fermi energy E F ≈ 100 meV and free carrier density of N c ≈ 10 12 cm −2 (see Equation (4)), [59,63,85] available for the interaction with the THz fields via intraband conductivity, and leading to a typical THz power absorbance on the order of A intra ≈ 0.1 at room temperature (see, e.g., refs. [59,63]). Second, as explained in Section 3, the Pauli-blocking of interband THz transitions at elevated temperatures strongly suppresses the interband absorption of graphene, leading to a negligible interaction with the THz fields. The room-temperature power absorbance of graphene at 1 THz frequency is only A inter = πα(f FD, conduction − f FD, valence )| 1THz, 300K ≈ 10 −3 , which is ≈100 times smaller than the typical intraband THz absorbance for doped graphene of A intra ≈ 0.1. Here, α is the fine-structure constant (see Section 3). Therefore, THz nonlinearities related to interband transitions in graphene can only be observed in (almost perfectly) intrinsic graphene, and only at cryogenic temperatures. Indeed, nonlinear multi-photon-like induced THz absorption in near-intrinsic graphene was recently demonstrated under such conditions, originating from a combination of coherent and incoherent inter-and intraband excitation mechanisms. [62] However, this experiment required a multilayer sample (45 layers) in order to significantly enhance the light-matter interaction strength.
At the same time, the intraband transitions in graphene remain active even at room temperature and, as mentioned above, lead to strong THz absorbance with A intra ≈ 0.1. As a result, the interaction of graphene with THz waves via the intraband conductivity mechanism remains much stronger in both linear and nonlinear regimes, as compared to the interband THz transitions. [9,13,55,59,63,111] Figure 4 shows the state-of-the-art of the experimentally measured third-order nonlinear susceptibility χ (3) of graphene, in the spectral range from visible to THz. A significant increase of χ (3) with the fundamental wavelength of the driving light field is observed, yielding the largest nonlinear coefficient of χ (3) ≈ 10 −9 m 2 V −2 at the THz frequencies (case [g] in Figure 4). This value of χ (3) of graphene is several orders of magnitude larger than that of other materials in any spectral range, possibly making graphene the most nonlinear material known to date. We note that even when normalized to material thickness (in case of graphene d = 0.3 nm) [114] and using the sheet values for nonlinear coefficients χ (3) 2D = χ (3) d, the nonlinearity of graphene still by far exceeds that of other known highly nonlinear systems: χ ≈ − 10 2D (3) 23 m 3 V −2 for highly nonlinear quantum well structures with d of a few nm. [117,118]

Experimental Observation of THz Nonlinearities in Graphene
In the following, we will highlight the exemplary experimental results and measurement approaches that led to the state-ofthe-art understanding of the THz nonlinearity in graphene.

Experimental Techniques
Various techniques and different sources of intense THz fields have been employed for nonlinear THz spectroscopy of graphene. The most commonly used strong-field THz sources are table-top sources based on optical rectification of amplified laser pulses in various nonlinear crystals, [55,[57][58][59]62,119,120] or based on air-plasma THz generation method. [121] Most of such experiments represent the nonlinear (variable field strength) version of THz time-domain spectroscopy (THz-TDS), THz pump-THz probe, or THz pump-optical probe experiments, and are performed with spectrally broadband, near single-cycle THz pulses. A comprehensive description of high-field tabletop THz sources and nonlinear THz spectroscopic technics can be found in refs. [61,122].
Further, an important role in understanding of the IR and THz optical properties of graphene was played by THz sources at large-scale accelerator-based facilities, such as IR/THz free electron lasers (FELs) [123] and superradiant THz undulator facilities, [124] providing spectrally dense, narrowband, and broadly tunable output. The most common experiment at FELs has been based on degenerate one-color pump-probe spectroscopy. In a series of careful experimental campaigns, the induced transient differential transmission change in graphene has been measured on Fourier-limited few picosecond timescales over a wide range of frequencies stretching from few THz to midinfrared frequencies. [9,114,125,126] These experiments revealed important insights into the carrier dynamics in graphene, including the pioneering observation of inter-and intraband transitions in doped graphene. [9] Most recently, a prototype superradiant THz undulator facility TELBE [124] was used for the pioneering demonstration of THz HHG in single-layer graphene. [63] The key for the success of this experiment was a unique combination of emission and detection characteristics of TELBE: (i) spectrally pure, multicycle, phase-stable strong-field pulses of THz and sub-THz radiation delivered at exceptionally high repetition rates, and (ii) novel THz detection schemes based on precise femtosecond timing of the superradiant undulator facility to external laser sources. The former allowed for exceptionally good statistics on the measurement, and the latter enabled a particularly sensitive probing of nonlinear THz signals in the time-domain with subcycle temporal resolution and a dynamic range as high as 10 6 .

Key Observations of Nonlinear THz Absorption in Graphene
Here, we will present the key observations of nonlinear THz absorption in graphene, as well as several results of optical pump-THz probe spectroscopy that have been crucial for the understanding of the nature of THz nonlinearity of graphene.
In 2011, Hwang et al. [54,55] reported the pioneering observation of the THz nonlinearity of CVD-grown graphene, using intense single-cycle THz pulses with a maximum fluence of 190 μJ cm −2 corresponding to a peak electric field strength of ≈100 kV cm −1 , see Figure 5a-d. As experimental methods in this study, nonlinear THz-TDS and THz pump-THz probe spectroscopy were used, as depicted schematically in Figure 5a. In this work, a THz-field-induced transparency of doped graphene, that is, an increase of THz transmission through graphene with increasing the THz pump field, shown in Figure 5b,c, was demonstrated for the first time. The nonlinear transmission enhancement was attributed to saturable absorption effects related to a decrease in the intraband conductivity of the doped graphene sample. The measured power-dependent THz transmission of graphene was described, as shown in Figure 5c [112] [b] Kumar et al., [43] [c] Soavi et al., [46] [d] Hendry et al., [48] [e] Kundys et al., [113] [f ] König-Otto et al. (Landau-quantized graphene in a magnetic field), [114] and [g] Hafez et al. [63] Whereas in the NIR (left), the nonlinearity mechanism is typically associated with coherent electron dynamics, in the THz (right) the mechanism is noncoherent and thermal. For intermediate wavelengths, different mechanisms have been proposed, such as plasmon nonlocal effects [e], and thermally induced plasmon shifts [115,116] (no susceptibility quoted).
a phenomenological saturable transmission function yielding a saturation fluence of 20 μJ cm −2 . An attempt to reproduce the data quasi-quantitatively using a Drude model of the THz intraband conductivity of graphene (see Equation (7)) required variation of the electron momentum scattering time as a free fitting parameter: The trend of reduction of THz conductivity of graphene could be reproduced by reducing the electron scattering time in graphene with increasing the driving THz field.
Adv. Optical Mater. 2020, 8,1900771 Figure 5. THz field-induced transparency in monolayer graphene: a) Schematic of varying high-field THz-TDS and intense-THz-pump/weak-THzprobe spectroscopy of graphene on a substrate; techniques followed by Hwang et al. [55] b) Spectral field transmission increasing with the THz fluence, with the temporal THz fields (left inset) and the corresponding spectra (right inset) transmitted through the graphene sample and the bare substrate, c) The power transmission versus the THz pump fluence; symbols for the experiment and solid red line for the saturable absorption fit, d) Transient differential transmission of the THz probe beam transmitted through THz-pumped graphene as a function of the pump-probe delay time, showing ultrafast (less than 100 fs) transient increase in transmission followed by a subsequent biexponential relaxation over a picosecond time range (b-d: Adapted with permission. [55] Copyright 2013, American Chemical Society), and e) suppression of graphene conductivity with THz field for all the THz frequencies in test demonstrating the observed transparency (Reproduced under the terms of the Creative Commons CC-BY 4.0 License. [59] Copyright 2015, Springer Nature).
In the THz pump-THz probe experiment presented in the same work, [54,55] the dynamics of THz absorption bleaching in graphene was demonstrated, with quasi-instantaneous increase in THz transmission followed by a biexponential relaxation on ≈ 200 fs and ≈picosecond timescales (see Figure 5d). The pioneering observation of THz nonlinearity of doped graphene by Hwang et al. [54,55] was followed by further experiments on nonlinear THz spectroscopy of graphene [56,58,59,62,119,120] using table-top intense near single-cycle THz pulse sources, and similar observations of THz-field-induced transparency (THz absorption bleaching) in both CVD [56,58,59,111] and epitaxial [57,119] graphene were reported. The physical mechanism of THz absorption bleaching was mainly attributed in all these works to nonlinear intraband conductivity modulation (suppression) by the intense THz driving field (see Figure 5e). However, the precise mechanism of THz conductivity suppression remained under debate at the time.
In 2014, Bowlan et al. [62] observed induced absorption and nonlinear spectral modification of near single-cycle broadband strong-field THz pulses in near-intrinsic epitaxial 45-layer graphene sample at a cryogenic temperature, which was assigned to a complex process initiated by nonlinear THz interband transitions in graphene, vanishing with increasing the sample temperature. In a more recent work, four-wave mixing in multilayer (50 layers) epitaxial graphene in strong magnetic field, and driven by intense far-IR field at 19 THz was reported by König-Otto et al. [114] In 2015, Mics et al. [59] published a combined experimentaltheoretical study presenting a thermodynamic picture of the THz nonlinearity in graphene. This work largely reconciled the previous observations of THz nonlinearities in graphene, and provided a quantitative description of the mechanism of the strong modulation of intraband conductivity in the THzexcited graphene using common conservation laws. The basic mechanism of THz nonlinearity of graphene, which will be described in detail in Section 4.3, is as follows. The absorption of the THz field by free electrons of graphene leads to quasiinstantaneous increase in the electron temperature, accompanied by a downshift in the chemical potential, which in turn suppresses the intraband conductivity of graphene. [58,62,81] This rather simple thermodynamic model of the THz nonlinearity of graphene by Mics et al. [59] is currently believed to be the most accurate one, showing very good agreement with the experimental results on picosecond-timescale THz nonlinear optics of graphene.
With one-color pump-probe spectroscopy using a tunable FEL, in 2011 Winnerl et al. [9,126] studied ultrafast carrier dynamics in quasi-intrinsic (lightly doped) multilayer (70 layers) epitaxial graphene grown on SiC at a range of photon energies from 10 meV (≈2.5 THz) to 250 meV (≈62.5 THz), and made the pioneering observation of, and distinction between, the inter-and intraband transitions in graphene. For photon energies larger than twice the Fermi energy (i.e., ℏω > 2E F ), pump-induced occupation of energy states that are ω 1 2 above the Dirac point via interband transition leads to an absorption bleaching of the probe beam similar to semiconductors, resulting in a positive differential transmission signal (increase in transmission), as shown in Figure 6a-c for the case of ℏω = 30 meV. Increasing the pump fluence, as expected, results in a larger differential transmission ΔT/T signal. We note here the k-space distribution of the photoexcited carriers in graphene is initially anisotropic. This anisotropy was revealed in pump-probe FEL experiments by examining the polarization relation between the pump and probe beams, [93,126,127] as shown in Figure 6d,e. It is also observable for a broad range of IR excitation photon energies and even at THz frequencies, [128] for as long as a sufficient pump fluence is available, for both inter-or intraband excitations. This anisotropy vanishes via carrier relaxation by phonon emission within sub-picoseconds or on a longer timescale, depending on the energy of the excited state relative to the optical phonon energies in graphene. [93,126] Excitations in lightly doped graphene below the optical phonon energy (≈200 meV) are accompanied by a longer preservation of such an anisotropy over picoseconds due to quenching of the optical phonon emission. [126,127] In contrast to the case of high photon energy inducing interband transitions, pumping at lower photon energies (i.e., ℏω < 2E F ) in one-color pump-probe FEL experiment leads to interband absorption of the probe beam: negative probe differential transmission ΔT/T 0 for ℏω = 20 meV in Figure 6a. This is a result of smearing out the carrier distribution in the band by heating the carriers via free-carrier absorption (intraband dynamics) of the pump beam, thus opening channels for the probing photons of the same energy to induce interband transitions, as illustrated in Figure 6c. We note that this latter observation of stronger THz probe absorption in excited graphene seems to be in qualitative disagreement with the observations of THz probe transmission enhancement by Hwang et al. [55] shown in Figure 5d. This discrepancy can nevertheless be accounted for by the much lower doping concentration of the 70-layer graphene sample, and the relatively large photon energies ℏω = 20 meV = 4.9 THz employed in the study by Winnerl et al. Indeed, according to the explanation in Section 3, for such high photon energies the Pauli-blocking of interband transitions at room temperature becomes relaxed, thus permitting measurable interband THz absorption in a low-doped multilayer sample.
Let us now consider the optical pump-THz probe experiments, which produced a bulk of valuable knowledge, and significantly contributed to the present understanding of the nature of THz nonlinearity of graphene. In the following, we will notice the similarities between the THz response of graphene to strong THz excitation and to optical excitation.
Optical excitation of graphene significantly modifies the electronic occupation of its band structure. Absorption of an optical pulse by graphene leads to the creation of initial electron-hole pairs, which thermalize among all charge carriers within ≈100 fs, leading to a transient electron distribution with an elevated temperature, that is, to carrier heating. This ultrafast electron heating in photoexcited doped graphene, that is, the smearing out of electron Fermi-Dirac distribution in the energy-momentum space, was first directly revealed in 2013 in a pioneering time-resolved angular-resolved photoemission spectroscopy (tr-ARPES) experiment by Gierz et al. [100] performed with both inter-and intraband IR excitations.
An important excursion must be made here. If the electronic temperature becomes high enough (≈3000 K), this photoinduced carrier heating also leads to a weak (due to impedance mismatch with the far-field), very short-living (<100 fs), yet observable, spectrally broad optical emission, essentially representing the blackbody radiation by the photogenerated hot electron population in graphene. This effect was discovered in 2010 by Lui et al. [130] in pristine graphene under 800 nm optical excitation. In particular, the high-frequency tail of the blackbody emission, with photon energies exceeding the excitation photon energy of 1.5 eV, was observed with sensitive photodetectors in the range 1.7-3.5 eV. [130] This optical emission due to the blackbody radiation is a collective process of the hot carrier electron distribution, and is thus very different from the conventional semiconductor-type radiative electron-hole recombination. For example, this blackbody optical emission, unlike electron-hole recombination in semiconductors, does not lead to a reduction in the density of excited carriers in graphene. The presence of such blackbody emission in photoexcited graphene, [130] and the previously discussed time-resolved ARPES measurements showing the ultrafast increase in electronic temperature in photoexcited graphene, [100] together provided convincing proof of the occurrence of photoinduced carrier heating in graphene.
For near-intrinsic or lightly doped graphene with E F of the order of a few tens of meV, [84,119,120] ultrafast photoexcitation leads to an increase of the THz-probed graphene conductivity as a result of the increase in free carrier density, and observed as negative differential transmission of the probe THz pulse. [131,132] This effect is similar to above-bandgap photoexcitation of semiconductors.
However, if the photoexcited graphene is significantly doped, the situation becomes quite different. For a typical doping level with E F > 100 meV, several groups [14,15,84,119,120,129,131,132] observed a reduction of THz probe absorption of doped graphene as a result of photoexcitation. In Figure 7, we show the exemplary results by Tielrooij et al. [129] and Jensen et al. [84] This negative photoconductivity is in a way counterintuitive, as it shows that ultrafast optical interband excitation in fact makes doped graphene more THz-resistive. In a combined experimental-theoretical study in 2013, Tielrooij et al. [129] assigned this phenomenon to "hot-carrier multiplication," a cascaded electron-electron scattering process, triggered by the energetic photoexcited electrons, which scatter with the background electrons in the Fermi sea and increase their energy. The microscopic theory for such hot-carrier multiplication by multiple electron scattering in graphene was presented by Song and Levitov in the work of Song et al., [98] where it was also demonstrated that the rate of such electron-electron scattering generally largely exceeds the electron-phonon scattering rate. Such a photoinduced hot-carrier multiplication was connected in the work of Tielrooij et al. [129] (see its Supporting Information) to very fast electron thermalization and hence to the increase of the electronic temperature in graphene, leading to the reduction of its chemical potential and conductivity. These studies, explaining the negative photoconductivity of doped graphene, became quite important also for the understanding of Adv. Optical Mater. 2020, 8,1900771 Figure 6. One-color pump-probe spectroscopy of epitaxial multilayer graphene using narrow-band tunable free-electron laser pulses: a) Differential transmission of the probe beam as a function of the pump-probe delay time at various pump fluences for two cases of photon energy; ℏω = 30 meV > 2E F and ℏω = 20 meV < 2E F (Reproduced with permission. [126] Copyright 2017, Wiley-VCH, while the original content of the figure is published in Ref. [9]), b) model calculations of the corresponding real conductivity at different electron temperatures over a range of photon energies, showing pump-induced absorption of the probe beam when ℏω < 2E F and induced transparency when ℏω > 2E F (Reproduced with permission. [126] Copyright 2017, Wiley-VCH), c) schematic explanation of the revealed responses in (a) and (b) (Reproduced with permission. [9] Copyright 2011, American Physical Society), d) schematic of anisotropic distribution of photoexcited carriers accumulated in the direction perpendicular to the polarization direction of the pump beam (Adapted with permission. [127] Copyright 2014, American Chemical Society), e) differential transmission signal revealing the anisotropy described in (d) using one-color pump-probe spectroscopy in which parallel and perpendicular polarizations of the pump and probe beams were compared (Reproduced with permission. [93] Copyright 2016, American Physical Society).
THz-induced nonlinearities in graphene. In Section 4.3, we will show that the final result of both intense THz (with an electric field strength exceeding 10 kV cm −1 , fluence of ≈2 μJ cm −2 ) and optical excitation (usually with a fluence exceeding 1 μJ cm −2 ) of doped graphene is, in fact, equivalent-an increase of the temperature of its background electron population, which in turn reduces its intraband, that is, THz conductivity. [59,63,84,129] The doping concentration in graphene was found to play a significant role in defining the THz response of graphene to both optical and intense THz excitations. Razavipour et al. [111] reported the effects of the carrier concentration in graphene on the THz field-induced nonlinearity (only absorption bleaching), using gated graphene samples in which the doping concentration was controlled via electrochemical gating, as depicted schematically in Figure 8a. It was found that the THz field-induced transparency in graphene (or the THz nonlinearity in general) becomes more pronounced with increasing the doping concentration, as shown in Figure 8b,c. This indicates that the nonlinear THz field effects in graphene are dominated by intraband conductivity mechanisms, so that the higher the doping concentration the larger the possible modulation of the intraband conductivity by the THz field is.
When the free carrier concentration in the sample was controlled by the combination of electrical doping by gating (static doping) and ultrafast optical excitation, the situation became more complex. [84,119,120,131,132] Depending on the initial doping density and the THz probe field strength, the THz photoconductivity was found to be either negative (THz absorption bleaching) or positive (enhancement of THz absorption), as shown in Figure 8d-f. In lightly doped graphene with E F on the order of a few meV, increased THz absorption in photoexcited graphene was noticed when probing by weak THz field (see Figure 8e), similar to the results in ref. [84] shown in Figure 7c, indicating the semiconductor-like THz response of near-intrinsic graphene. With increase in the THz probe field strength, THz absorption bleaching with a crossover to negative THz photoconductivity (positive THz differential transmission) was observed. On the other hand, for strongly doped graphene the THz photoconductivity always remained negative, indicative of THz absorption bleaching (positive THz differential transmission) as a result of combined action of THz and optical excitation, as shown in Figure 8f.
In a more recent work, Tomadin et al. [133] significantly extended the understanding of the effect of photoexcitation on graphene based on a combination of microscopic modeling and optical pump-THz probe spectroscopy. Interestingly, it was found that depending on the Fermi energy of the sample, the ultrafast electronic thermalization in photoexcited graphene can lead Adv. Optical Mater. 2020, 8,1900771 [129] Copyright 2013, Springer Nature), b) THz differential transmission of the peak of the THz probe field as a function of the pump-probe delay time for two cases of optical-pump photon energies maintaining the same photoexcited carrier density of ≈10 11 cm −2 , showing an increase in the THz transmission of graphene after photoexcitation, corresponding to a suppression of the photoconductivity (Adapted with permission. [129] Copyright 2013, Springer Nature). c) Crossover from negative to positive photoconductivity in graphene with decreasing free carrier concentration. Negative photoconductivity in doped graphene and positive photoconductivity in near-intrinsic graphene at fixed photoexcitation level. A sample with a back-gate was used in order to control the position of the Fermi level in graphene (Adapted with permission. [84] Copyright 2014, American Chemical Society). not only to "hot-carrier multiplication," that is, electron heating with the total number of free carriers conserved (case of larger E F > 0.1 eV), but under certain conditions also to a "real" free carrier multiplication via interaction with the valence band states (case of smaller E F < 0.1 eV). Carrier multiplication in graphene was first predicted by Winzer et al. [134] and observed experimentally by Plötzing et al. [135] Ref. [133] furthermore provides insights into the mechanism leading to the reduction of the graphene conductivity with increasing electron temperature, based purely on electron kinetics effects, including the effect of screening, and considering that momentum scattering is dominated by longrange Coulomb scattering with impurities outside graphene.
An interesting observation was made by Tani et al. [121] in a THz pump-optical probe experiment, performed using ≈30 fs single-cycle THz pump pulse and NIR probe of 800 nm wavelength (1.55 eV photon energy). A standard doped CVD-grown graphene was used in this work. It was shown that the excitation of graphene with the THz pulses with electric field of up to 300 kV cm −1 resulted in an ultrafast (≈150 fs) transient increase of the optical probe transmission through graphene, as shown in Figure 9. This effect was assigned to absorption bleaching of the optical probe by Pauli-blocking of the NIR interband transitions at 1.55 eV photon energy under intense ultrafast THz excitation of graphene. The authors attributed the electron occupation of states around 0.775 eV above the Dirac point to impact-ionization effects. The impact ionization is the process of electron-hole pair generation resulting from the loss of kinetic energy of another highly energetic electron or hole, which was assumed to be ballistically accelerated by the strong pump THz field. It was suggested that this process leads to carrier multiplication, allowing for the occupation of higher energy states around 0.775 eV above the Dirac point, in turn bleaching the interband optical absorption of the 1.5 eV probe pulse in graphene. The impact ionization in strong (up to 300 kV cm −1 ), ultrafast (≈30 fs) THz field pulses as observed by Tani et al. [121] increases the free carrier density in graphene. This is an additional effect to the THz-induced electron heating-the main effect of interaction of doped graphene with intense THz fields, which as such conserves the free carrier density (see Section 4.3 for details). To the best of our knowledge, the impact ionization leading to carrier multiplication has not been observed in the experiments dealing with "slower" (>1 ps long) strong-field THz pulses as used in the experiments of refs. [56][57][58][59]63,111,119].
Besides standard "continuous" graphene samples, nonlinear THz spectroscopy was also performed on patterned graphene, featuring the confinement of THz field, and thus leading to enhancement of the THz nonlinearity. For example, it was theoretically predicted [136][137][138] and experimentally verified [51,115,116] that enhancement of the nonlinear interaction of the light field in graphene can be achieved by local field confinement via resonant plasmonic absorption if the graphene is patterned in subwavelength structures such as ribbons. Using THz FEL, Jadidi et al. [115,116] demonstrated nonlinear THz absorption at plasmonic resonance in lithographically patterned epitaxial graphene ribbons (see Figure 10a), following one-color pump-probe scheme with THz field polarization perpendicular to the patterning direction. The plasmonic resonance frequency (see Figure 10b) was defined by the dimensions of the graphene pattern, as well as by the electron density and temperature. The THz pump at resonance resulted in two consequent effects: (i) a rise in the electron temperature by the THz excitation enhanced due to the THz field confinement, and consequently (ii) a redshift in the resonance absorption. The probe experienced absorption Adv. Optical Mater. 2020, 8,1900771 [111] Copyright 2015, American Physical Society). d) Schematic of optical pump-induced dynamics probed by intense THz field, and e,f) peak of THz differential transmission signal ΔT/T 0 | peak of photoexcited graphene as a function of the THz probe field for lightly doped graphene pumped by 29 μJ cm −2 optical pump fluence (Adapted with permission. [120] Copyright 2015, AIP Publishing), highly doped graphene at three optical pump fluences of 29, 79, and 137 μJ cm −2 (Adapted with permission. [119] Copyright 2015, American Physical Society). Figure 9. THz pump/optical probe of single layer graphene: a) Schematic demonstration of the THz field induced dynamics in hole-doped CVD graphene, showing ballistic acceleration imposing impact ionization from the valence band to the conduction band leading to occupation of states that block further interband transitions by the optical probe beam, b) the THz pump field with a peak field of 300 kV cm −1 with the corresponding spectrum shown in the inset, c) negative differential optical density corresponding to enhanced transmission of the optical probe beam as a function of the pump-probe delay time, d,e) experimental data (symbols) of the optical density change due to the intense THz pump as a function of the pump-probe delay time, and model calculations (solid lines) based on the Boltzmann equation with scattering mechanisms including impact ionization and Auger recombination for THz field strengths of 170 and 300 kV cm −1 , respectively, f,g) calculated carrier density following the impact ionization and Auger recombination, as a function of the pump-probe delay time, and h) the peak field dependence of the peak of the differential optical density, showing a scaling with THz 2.5

E
and a better fit with the model only when impact ionization and Auger recombination were considered (Reproduced with permission. [121] Copyright 2012, American Physical Society). . The normalized spectrum in red corresponds to one of the FEL pulses (here at resonance) employed in the degenerate pump-probe spectroscopy (with arrow indicating relevant axis) (Replotted with data provided by the authors and with permission. [116] Copyright 2019, American Chemical Society). c) THz differential transmission at resonance of the probe beam as a function of the pump-probe delay time for various THz pump fluences (Reproduced with permission. [116] Copyright 2019, American Chemical Society), d) The peak differential transmission signal as a function of the pump fluence; colored symbols for the experiment and a solid black line for the F fit (Reproduced with permission. [116] Copyright 2019, American Chemical Society), e) the differential transmission signal at various frequencies around the resonance frequency showing a sign reversal when the probe frequency is tuned to be below resonance (Reproduced with permission. [116] Copyright 2019, American Chemical Society), and f) the peak of the differential signal in (e) as a function of frequency (Reproduced with permission. [116] Copyright 2019, American Chemical Society). bleaching, with positive differential transmission signal ΔT/T 0 increasing as the square root of the pump fluence, as shown in Figure 10c,d. The relatively small pump fluence employed here, as compared to that employed in the studies with continuous planar graphene, [9,55,59] showed a significant enhancement in the THz field-induced nonlinearity by resonant plasmons in patterned graphene. On the other hand, the redshift of the plasmon resonance due to the rise in the electron temperature induced by the THz pump led to increased absorption of the probe beam when probing at frequencies lower than the original resonance frequency (negative ΔT/T 0 in Figure 10e,f).

The Thermodynamic Model of Nonlinear Intraband THz Response of Graphene
The observations of the nonlinear THz response of graphene, described in Section 4.2.2, can be summarized as follows. For doped graphene, excitation with intense THz fields (field strength above ≈10 kV cm −1 , fluence of the order of ≈μJ cm −2 ) leads to a reduction of the THz conductivity of graphene. The effect of THz excitation of doped graphene is thus the same as the effect of optical excitation (also summarized in Section 4.2.2), which similarly reduces the THz, that is, intraband, conductivity of graphene.
In the following, we will describe the basic mechanism of such a THz-driven conductivity suppression, following the thermodynamic model of ultrafast charge transport in graphene presented by . [59] This model, resting on basic conservation laws, self-consistently treats the THz-induced electron heating in doped graphene during its interaction with the incident THz field, dynamically modifying the propagation conditions for the THz field itself. The key physical processes described by the thermodynamic model of THz nonlinearity of graphene are summarized in Figure 11.
We note here that most equations of this model are only strictly valid in the Sommerfeld expansion regime, that is, at electron temperature being below the Fermi temperature: T e < T F , where T F = E F /k B . For highly doped graphene on SiO 2 with typical Fermi energy of 0.25-0.4 eV, the Fermi temperature is T F = 3000-4800 K. For practical purposes, however, the approximation T e < T F in the description of nonlinear interaction of Adv. Optical Mater. 2020, 8, 1900771   Figure 11. The thermodynamic model of THz-driven electron transport in graphene: a) Schematic for the THz field-induced electron dynamics via intraband THz conductivity (THz absorption by the Fermi surface electrons) showing intraband heating enabled by e-e scattering that allows electrons in the band to thermalize among themselves, which is accompanied by a reduction of the chemical potential to keep the carrier density constant (Reproduced under the terms of the Creative Commons CC-BY 4.0 License. [59] Copyright 2015, Springer Nature), b) the population density before (blue) and after (red) the THz excitation (Reproduced with permission. [63] Copyright 2018, Springer Nature), c) the corresponding Fermi-Dirac distribution before and after the THz excitation showing a smearing-out of the hot carrier distribution in the band accompanied by a downshift in the chemical potential to maintain the total carrier density (the area under the curves in (b)) constant, and d) illustration of spectral weight conservation for inter-and intraband transitions in graphene. Reduction of chemical potential for hotter electron population leads to increase of the spectral weight for interband transitions (as a result of Pauli-unblocking), which is compensated by the reduction of the spectral weight of intraband conductivity in graphene (c,d: Reproduced with permission. [82] Copyright 2015, IOP Publishing).
arbitrary THz field with graphene remains reasonably accurate. Most of the nonlinear response in graphene is in fact achieved in the regime of T e < T F , and the effect of further increase in electron temperature appears smaller in comparison. For higher electron temperatures, leading to chemical potential approaching the Dirac point, the accuracy of the thermodynamic model by Mics et al. [59] can be improved by including the valence band states, as shown by Frenzel et al. [132] Let us first consider the energy conservation in the interaction of THz radiation with doped graphene. The application of the THz field E THz to graphene leads to generation of THz currents according to Ohm's law, j = σ intra E THz , where σ intra is the intraband conductivity of graphene, described by Equation (6). Generally, the conductivity is directly proportional to the power absorption coefficient of the material α = σ′/cnε 0 , where σ′ is the real part of conductivity (see Equation (6)), n is the real part of refractive index, c is speed of light in vacuum, and ε 0 is vacuum permittivity. The excitation of THz current j in doped graphene thus facilitates the transfer of energy from the absorbed THz field to the free conductive electrons in graphene. In the approximation of sharp Fermi edge, which is convenient for this discussion, one can assume that this energy is initially confined to the electrically conductive states residing at the Fermi level E F . The linear energy-momentum dispersion in graphene given by Equation (2), apart from the high electron mobility due to fundamental absence of Coulomb binding, yields yet another unique property: perfect energy and momentum conservation for collinear electron-electron scattering. In a simple intuitive picture, given available phase space, any two Dirac electrons in graphene can scatter by conserving their total energy and momentum: One electron increases its energy and momentum by ΔE and Δk, respectively, while the second electron decreases its energy and momentum by the same amounts (see Figure 11a). Detailed considerations indeed showed that collinear scattering between two graphene electrons is the most probable, and therefore the fastest electron scattering process. [98,99,139] This is in contrast to materials with quasi-parabolic energy-momentum dispersion, such as semiconductors, where typically more than just two participating particles are needed in order to fully conserve both energy and momentum in a scattering event, thus reducing its probability and consequently slowing the scattering rate down.
It is via these collinear electron-electron interactions in graphene, that the energy of the absorbed THz field, initially confined to the conductive states at the Fermi edge, is being redistributed across the entire carrier population of graphene, becoming its collective kinetic energy-the electronic heat. The process of such an energy redistribution, the electron thermalization, is known to happen on an ultrafast timescale of <100 fs for both intra-and interband excitation of doped graphene. [100,[140][141][142][143] This thermalization is much faster than the oscillation cycle of THz waves, which is on the order of picoseconds. As a result, one can treat the thermalization dynamics as an instantaneous process on the timescale of THz-graphene interaction, which leads to the following thesis: (1) Energy conservation: The THz energy absorbed by graphene electrons via intraband conductivity is instantaneously converted into electronic heat with perfect conversion efficiency, raising the electronic temperature.
Since the electrons are assumed as internally thermalized at all times, and hence are described at all times by Fermi-Dirac  with a certain temperature T e and chemical potential μ, the free electron conductivity of a THz-driven graphene at all times can be described by Equation (6) depending on the above parameters. Quantitatively, the total thermal energy density (energy per unit area) of the graphene electronic system excited by THz field can be expressed as ∫ ε ε ε µ ε is the density of states. Here, Q 0 is the electronic thermal energy before the THz field excitation, that is, with the chemical potential μ equals to the initial Fermi energy at the initial temperature (e.g., room temperature). With the THz field absorption, an additional heat ΔQ is added, thus leading to a rise in T e , as illustrated in Figure 11a.
Second, we need to impose a condition on electron density conservation for the intraband THz excitation of graphene: at any electron temperature This requires that the rise in temperature T e due to the THz deposited heat ΔQ must be accompanied by a reduction in the chemical potential μ, as shown in Figure 11b,c (see also ref. [82]). The condition of carrier density conservation thus leads to the second thesis of the thermodynamic model: (2) Carrier density conservation: Increase of electron temperature leads to reduction of chemical potential in THz-excited graphene.
Now that we established that the intraband THz absorption in graphene leads to an increase of electron temperature accompanied by a downshift in the chemical potential, we arrive at a qualitative understanding of the reduction of THz conductivity in THz-excited graphene, based on the spectral weight conservation rule, illustrated in Figure 11d. This rule was first applied to the treatment of THz photoconductivity in graphene by Frenzel et al. [132] Indeed, the increase of the electron temperature, resulting in the downshift of the chemical potential μ, leads to unblocking of spectral density for previously Pauli-blocked interband transitions at the (infrared) photon energy of ħω > |2μ|. This increases the spectral weight for interband transitions for graphene with hotter electron population, whereas the total spectral weight of both interband and intraband transitions must remain conserved. This brings us to the third thesis of the thermodynamic model.
(3) Spectral weight conservation: The downshift of chemical potential in graphene with higher electron temperature leads to an increase of the spectral weight for interband transitions, compensated by a decrease of the spectral weight for intraband transitions. The latter reduces the THz conductivity in graphene with higher electron temperature. Now we consider the dynamics of the processes underlying the thermodynamic model of conduction in graphene. While the energy transfer from the driving THz field to the free electron population of graphene, raising the electronic temperature, occurs quasi-instantaneously on sub-100 fs timescale, the consecutive electron cooling occurs on a slower, few-picosecond timescale of phonon emission, as schematically depicted in Figure 11b. This dynamics of electron cooling after intraband excitation of doped graphene was measured in the pioneering tr-ARPES experiments by Gierz et al., [100] see Figure 12a. The electron cooling dynamics features three distinct timescales ranging from ≈10 fs to 1 ps, which are assigned to the electron relaxation via optical and acoustic phonons, respectively (see also Table 1). The disparity between electron heating (sub-100 fs) and total cooling (several ps) timescales naturally leads to the heat accumulation in the electronic population of THz-driven graphene, ultimately resulting in the reduction of THz conduc-tivity in graphene on the timescale of the oscillation cycle of the THz field. Note how this dynamics of electron cooling, directly observed in tr-ARPES measurements, [100] roughly corresponds to the dynamics of THz conductivity recovery in a THz pump-THz probe experiments by Hwang et al. [54,55] (see Figure 5d).
Further, the tr-ARPES data by Gierz et al. [100] demonstrated saturation of the electron temperature in graphene with increasing intraband excitation fluence, see Figure 12b. Such Adv. Optical Mater. 2020, 8,1900771 [100] Copyright 2013, Springer Nature). c) Temporal evolution of electronic temperature in graphene with E F = 70 meV under quasi-single-cycle THz excitation of variable field strength (gray area indicates the modulus of driving THz field). d) Measured (symbols) and calculated (solid lines) nonlinear conductivity spectra of doped graphene, as the peak THz field increases from 2.3 kV cm −1 (yellow) to 120 kV cm −1 (black). e) Measured and calculated nonlinear absorption of doped graphene with E F = 70 meV, and the accumulated thermal energy density retained in the electronic system of graphene on the timescale of the interaction with quasi-single-cycle THz pulse, as well as the peak and average (in brackets) values of electronic temperatures as functions of THz peak field strength (arrows indicate relevant axis) (c-e: Reproduced under the terms of the Creative Commons CC-BY 4.0 License. [59] Copyright 2015, Springer Nature). a saturation indicates a fundamental nonlinearity in the intraband excitation of graphene: raising the electron temperature suppresses the intraband conductivity, thus limiting the energy transfer rate for the following portion of the excitation signal. This saturation occurs on top of the sublinear relation between the electron temperature and the incident power, which is caused by the electron heat capacity increasing with electron temperature. Note again the similarity in the saturation of tr-ARPES-measured electronic temperature in the work by Gierz et al. [100] (Figure 12b) and the saturable THz absorption in THz doped graphene from the pioneering work of Hwang et al. [54,55] (Figure 5c). Now that both the mechanism of energy transfer from THz field to graphene-the intraband conductivity according to Equation (6), and its dynamics-quasi-instantaneous electron heating and picosecond-timescale electron cooling (see Figures 11b and 12a) are established, one can numerically model the interaction of graphene with the driving THz fields in the time domain.
We note here that the intraband conductivity given by Equation (6) also contains an electron energy-dependent scattering time τ(ε). For graphene in full thermal equilibrium, that is, before strong THz excitation, this scattering time corresponds to the scattering time at the Fermi level τ(E F ), which can be determined experimentally in a linear THz-TDS experiment such as Figure 2f. Even though several electron scattering regimes are postulated for (hot) electrons in graphene, there is no full consensus on this topic in the literature until now. The two dominating momentum scattering mechanisms are known to be the long-range Coulomb scattering and short-range disorder scattering, for which the scattering times can be approximated as, respectively, τ (ε) = γε and τ (ε) = β/ε, where γ and β are linear proportionality constants (see, e.g., refs. [28,59] for discussion). Typically, long-range Coulomb scattering with τ (ε) = γ ε is believed to dominate for CVD-grown graphene deposited on a regular dielectric substrate (such as SiO 2 ), whereas this scattering channel can be suppressed in exfoliated graphene that is suspended or encapsulated by hexagonal boron nitride. Assuming that one of these scattering mechanisms dominates, the proportionality constant, whether γ or β, can thus be established experimentally from a linear THz-TDS experiment at room temperature as γ = τ(E F )/E F and β = τ(E F )E F , respectively. Both the single value for scattering time τ, and the Fermi energy E F can be readily determined in such a measurement (see Figure 2c-f and refs. [59,63,144]).
Once the dominating scattering mechanism for hot electrons, and hence the energy dependence of the scattering time, is assumed, one directly applies the conditions of energy and electron density conservation to the interaction of the driving THz field with graphene. The experimentally measured incident THz field is then numerically propagated in the time domain through the graphene, which at each point in time is characterized by instantaneous conductivity from Equation (6) dependent on instantaneous value of electron temperature T e and chemical potential μ, and with the energy-dependent scattering time τ(ε). This numerical propagation generates the transmitted THz field that can be analyzed in the same fashion as in nonlinear THz-TDS experiments, yielding the nonlinear, field-dependent conductivity spectra of graphene.
Such a numerical propagation can be performed fully numerically as is customary in traditional nonlinear optics (see, e.g., refs. [145,146] and references therein), or using convenient semi-analytical parameterization (see refs. [59,63] including the Supporting Information and Methods sections).
In Figure 12c,d, the results of such a numerical propagation of the experimentally measured incident THz pulse through the CVD graphene with E F = 70 meV and τ(E F ) = 140 fs are shown, assuming the long-range Coulomb scattering scenario for the hot electrons. [59] Figure 12c shows the calculated evolution of the electronic temperature in graphene, driven by near-single-cycle THz pulse of variable field strength, whose modulus is shown as a gray area. The evolution of electronic temperature T e generally reproduces the dynamics of the THz field. However, on the timescale of the light-matter interaction T e never reaches the initial temperature of 300 K, since the total cool-off of the electronic population is a few-picosecond timescale process, which is only complete after the main part of the driving THz field has gone (see ref. [100] and Figure 12a). This results in permanent electronic heat accumulation in graphene on the timescale of the interaction with the driving THz field, reducing the overall THz conductivity and absorbance of graphene, as the THz field is transmitted through it.
The numerical transmission of the THz field through graphene using the thermodynamic model yields the nonlinear THz spectra of graphene, shown in Figure 12d (solid lines) together with the experimental data (symbols), for the THz fields ranging from 2.3 to 120 kV cm −1 . [59] The conductivity in THz-driven graphene decreases with both the THz field strength and the frequency. It is evident that the agreement between the calculations and the entirety of the experimental data is quantitative, lending credence to the theses underlying the thermodynamic model of the nonlinear THz conduction of graphene. [59] Figure 12e shows the measured and the calculated power absorption of the graphene sample as function of the THz field strength, demonstrating a very strong saturable absorption effect, as has also been observed in many previous measurements (see, e.g., refs. [55,56,58,59,111] and Section 4.2 for discussion). The THz power absorption of doped graphene, measured here over the full timescale of THz-graphene interaction of several picoseconds, reduces from over A = 0.125 down to almost zero as the driving THz field increases from 2.3 to 120 kV cm −1 . This effect of very strong THz intraband absorption bleaching is in perfect agreement with the saturation of electron temperature in doped graphene under intraband excitation, observed in tr-ARPES experiments, [100] which, as mentioned before, in turn limits the energy transfer rate from the driving THz field to graphene electrons. In the same figure, the calculated accumulated electron heat, as well as the peak and average (in brackets) electronic temperature of graphene is presented. The energy balance shows that in this case about 15% of the absorbed THz energy remained within the electronic system of graphene on the timescale of interaction with the THz field, irrespective of the THz field strength. The calculated electronic temperatures in 1000s of Kelvin range, reached as a result of THz excitation, are in agreement with the measurements of electronic temperature in tr-ARPES experiments. [100] As we already mentioned in Section 4.2.2, the thermodynamic model of intraband electron dynamics in graphene is also applicable to the THz response of doped graphene under optical interband excitation. The tr-ARPES experiments with interband pump, [100,147,148] as well as optical pump-probe spectroscopy [139,149] showed a similar effect of increase of electron temperature in doped graphene subject to interband excitation. Indeed, also in the case of interband optical excitation the collinear electron-electron scattering of Dirac electrons will lead to at least partial redistribution of the energy of photoexcited electron-hole pairs to the background free carrier population, thus increasing its temperature and driving down THz conductivity of graphene. [84,129,133] As a result, THz conductivity of photoexcited graphene allows one to estimate the electronic temperature, and hence the efficiency of electron heating induced by optical excitation. We note here that the electron heating efficiency under optical interband excitation is not supposed to be perfect, unlike the THz-driven intraband heating efficiency. The reason is that the photoexcited electrons and holes usually possess enough energy in order to emit optical phonons of ≈200 meV energy, and hence to transfer their energy directly to the lattice of graphene, without exciting the background electron population. Optical phonon emission thus becomes the competing process to the background electron heating in graphene. In ref. [84], Jensen et al. used optical pump-THz probe spectroscopy to estimate the optically driven electron heating efficiency in photoexcited doped graphene. Using a numerical version of the model by Song and Levitov, ref. [129] and its Supporting Information, the photoconductivity spectra of photoexcited graphene were described, thus extracting the dependence of the electronic temperature and the total thermal energy of the electron gas in graphene on the fluence of photoexcitation. It was shown that in the case when the total energy of the unexcited electron gas ("depth of the Fermi sea") was significantly lower than the absorbed energy of the interband photoexcitation of graphene, all the energy of the photoexcited electrons and holes was converted into electron heat, without the direct excitation of the lattice via optical phonon emission. However, for higher excitation fluence and lower Fermi energy, the electron heating efficiency was found to decrease at the expense of the direct coupling of photoexcited electrons to the lattice. Tomadin et al. [133] further developed the findings of the work of ref. [84] using a more advanced microscopic model.

THz High-Harmonic Generation
As described in the previous section, the interaction of strong THz field E(t) with the free electrons in graphene results in strong suppression of graphene conductivity via carrier heating on sub-100 fs timescale, which is followed by the conductivity recovery via carrier cooling by phonon emission, occurring on the timescale of a few picoseconds (see Figure 11). As a result, during the interaction with the THz field, the instantaneous conductivity of graphene σ intra (E,t) becomes time-and field-strength dependent. This conductivity nonlinearity in turn leads to the generation of nonlinear current in graphene j(t) = σ intra (E,t)E(t) in response to the driving THz field E(t). If the driving THz field is monochromatic, that is, it oscillates at a defined frequency f, the induced nonlinear current in graphene will contain oscillations at higher order harmonics of the driving frequency. Due to the centrosymmetry of graphene, even-order harmonic generation is forbidden, and thus only odd-order higher harmonics 3f, 5f, 7f, … etc., will be generated. We note that the generation of high-order THz harmonics in graphene has long been proposed in theory (see, e.g., refs. [19,22,23]), yet the experimental demonstration of this long-sought effect has been achieved only recently.
In 2018, Hafez et al. [63] have demonstrated extremely efficient THz high-harmonics generation in graphene using the thermodynamic response of its free electrons to driving quasi-monochromatic THz fields, following the thermodynamic model by Mics et al. [59] The multicycle quasi-monochromatic THz fields with peak electric field up to 85 kV cm −1 at frequencies between 300 and 700 GHz, used as a pump signal, were generated by a superradiant THz undulator facility TELBE at Helmholtz-Zentrum Dresden-Rossendorf. [124] As a sample, a regular CVDgrown graphene deposited on a SiO 2 substrate was used, with the Fermi energy of E F = 170 meV, and the electron momentum scattering time at the Fermi level of τ(E F ) = 47 fs.
The quasi-monochromatic multicycle THz driving of graphene significantly enhances its effective nonlinear conductivity response, as compared to the single-cycle excitation demonstrated in refs. [55,56,58,59,111]. The key reason is that the electron cooling time, which is on the order of picoseconds, conveniently matches the oscillation period of a multicycle pump THz wave. This allows the graphene conductivity (and hence its absorbance) to efficiently recover in time for the next electric field cycle to arrive. Such an efficient recovery of conductivity of graphene, in turn, facilitates more efficient energy transfer from the THz pump wave to graphene electrons over the whole timescale of THz-graphene interaction, which for a multicycle THz pump wave can last 10s of picoseconds. Figure 13 shows the thermodynamic calculations, following the model by Mics et al. [59] of the THz-induced electron heat and the corresponding temporal modulation of graphene absorbance (and hence its conductivity), driven by the multicycle THz field oscillating at a fundamental frequency of 0.3 THz. This calculation particularly exemplifies the highly nonlinear nature of response of free electrons in graphene to THz excitation: the driving field strengths of 10 and 85 kV cm −1 lead to dramatically different effects on the electron temperature in graphene, and its absorbance. If the effect of 10 kV cm −1 excitation is minor, leading only to a slight perturbation of the equilibrium properties of graphene, the excitation with 85 kV/ cm −1 THz field strength results in a large THz-induced heat deposition into the system, and significant absorption and hence conductivity modulation featuring pronounced temporal oscillations observable twice per optical cycle of the driving field. This is a particular consequence of the centrosymmetry of graphene, where both the positive and the negative THz field oscillation half-cycle produces the same effect on its properties, in this case, suppressing its conductivity. Figure 14 shows the schematic of the experiment from the work of Hafez et al., [63] the spectrum of the driving THz field at the frequency f = 0.3 THz, as well as the spectrum of the THz field transmitted through graphene sample, with clearly visible higher odd-order THz harmonics up to the seventh order. The fundamental field at f = 0.3 THz is transmitted through the graphene with the transmission coefficient exceeding 0.9, which serves as a convenient reference for the conversion efficiency from the THz pump field to higher harmonics. These conversion efficiencies are remarkably high, especially for a material which is only one atomic layer thick: They demonstrate the values slightly less than 10 −2 , 10 −3 , and 10 −4 , for the third, fifth, and seventh THz harmonics, respectively. The total nonlinear field conversion efficiency therefore tends to ≈1%. The spectral range in this experiment was limited by the detector bandwidth, however, the signalto-noise ratio permitted by the TELBE THz source and the detection scheme [124] was not a limiting factor in this measurement.
In Figure 15a,b, the measured HHG signals in the time domain, as well as the corresponding results of the thermodynamic model calculations using the experimental pump wave are presented. The agreement between the measurement and the calculation is rather good: The thermodynamic calculation captures well all the key features of experimental signals, including their magnitude, distortion, and the shortening of the generated harmonic pulses with increase in their order. Figure 15c summarizes the experimental result of THz HHG in graphene by showing the electric field strength of the generated third, fifth, and seventh harmonics, as a function of the pump field strength (symbols). In the same figure, the results of the thermodynamic model calculation are shown, using the experimental pump field and basic parameters of graphene in full thermal equilibrium at 300 K as an input. The agreement between the data and the model is quantitative. From the observed dependencies in the regime of small-signal nonlinearities (pure power law, dashed lines), one can now estimate the effective nonlinear optical coefficients of graphene up to the seventh order: χ (3) eff ≈ 1.7 × 10 −9 m 2 V −2 , χ (5) eff ≈ 1.2 × 10 −22 m 4 V −4 , and χ (7) eff ≈ 1.74 × 10 −38 m 6 V −6 (see ref. [63] and its Methods section for details). These values are extremely large, and by many orders of magnitude exceed the nonlinear coefficients of any other known electronic material. This possibly makes graphene the most nonlinear material known to date.
Now that the thermodynamic model of THz response of graphene has demonstrated its high accuracy in reproducing Adv. Optical Mater. 2020, 8,1900771 Figure 13. THz field-induced heating and absorption modulation: a) Multicycle quasi-monochromatic THz driving field oscillating at 0.3 THz, b,c) the temporal THz induced electronic heat density (energy per unit area), and the THz power absorbance, respectively, of graphene driven by the intense multicycle THz pulse in (a) at two THz peak field levels of 10 kV cm −1 (the blue curves) and 85 kV cm −1 (the red curves), demonstrating the thermal response of the Dirac electrons of doped graphene to the intense THz driving field (Reproduced with permission. [63] Copyright 2018, Springer Nature).

Figure 14.
THz high-harmonic generation in single-layer graphene up to the seventh order at room temperature: a) schematic of the THz odd-order harmonics re-emitted from graphene driven by a multicycle THz field, b) experimental data of the THz field spectrum transmitted through the graphene sample (in blue) showing odd harmonics up to the seventh order relative to the spectrum of the driving field at a fundamental frequency of 0.3 THz (in red) (Reproduced with permission. [63] Copyright 2018, Springer Nature).
the THz nonlinearity in graphene to single-cycle [59] and multicycle [63] excitation, one can use it to predict the generation of THz harmonics of even higher orders. In Figure 16, the calculated spectra of response of graphene to the excitation at the fundamental frequency f = 0.68 THz are shown, predicting the appearance of harmonics up to the 13th order within the dynamic range of 10 6 , characteristic of the TELBE THz facility. The driving field amplitude in this calculation had a maximum of 1 MV cm −1 , which is a conservative estimate of the field strength required to enter the ballistic nonlinearity regime of coherent Bloch oscillations (see estimate provided in Section 4.1). We note here again, that the HHG output via the collective thermodynamic response of graphene electrons can be observed at the fields as weak as only 10s of kV cm −1 , which is 2-3 orders of magnitude weaker than that typically required for initiation of the coherent Bloch oscillation mechanism.

Conclusions and Future Perspectives
In this Progress Report, we have described the birth and evolution of the research area of THz nonlinear optics of graphene. We have shown that the strongest nonlinear response in graphene is facilitated via an intraband conductivity channel, and that the underlying physical mechanism is a very powerful thermodynamic, collective response of its free carrier population to the driving THz fields. We have summarized the fundamental aspects of a thermodynamic model of the THz nonlinearity of graphene, first described in the work by Mics et al., [59] which rests on basic conservation laws, and relies on perfect conversion of the THz energy absorbed via the intraband conductivity mechanism into electronic heat of the free electron population of graphene. This mechanism explains and quantitatively reproduces the experimentally observed strong THz saturable absorption and nonlinear wave conversion, including extremely efficient THz high harmonics generation in graphene at room temperature and under ambient conditions. As we show, graphene at THz frequencies possesses possibly the highest nonlinear coefficients of all materials we know to date.
We furthermore show that the thermodynamic mechanism underlying the enormous THz nonlinearity in graphene is active at rather moderate driving fields of only 10s of kV cm −1 . This THz field strength is an order of magnitude lower than the typical channel field of modern high-speed field-effect transistors, routinely operating at 100s of GHz frequencies. This is of high technological relevance since graphene is easily compatible with electronic technologies such as CMOS. [150,151] This creates the possibility to utilize CMOS-based sub-THz electronics as efficient and cost-effective pump sources for hybrid devices for Adv. Optical Mater. 2020, 8,1900771 Figure 15. THz high-harmonic generation in single-layer graphene up to the seventh order at room temperature-experiment versus calculations: a,b) the temporal THz fields of the fundamental driving field and the generated harmonics; experiment in (a) and calculations in (b). c) The dependence of the generated higher-harmonic field on the THz driving field; symbols for the experiment, solid black lines for the calculations, and dashed lines for the power-law fit of the low-signal nonlinearity regime below a driving field strength of 20 kV cm −1 (Reproduced with permission. [63] Copyright 2018, Springer Nature). Figure 16. Thermodynamic model calculations of THz high-harmonic generation up to the 13th order. Calculated spectrum showing high harmonics up to the 13th order within a dynamic range of 10 6 with respect to the transmitted field spectrum at the fundamental frequency (considering a fundamental frequency of 0.68 THz) (Reproduced with permission. [63] Copyright 2018, Springer Nature).