{'fraction_non_alphanumeric': 0.037615854860974164, 'fraction_numerical': 0.010648787221455335, 'mean_word_length': 4.342375559652357, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hot Electron Terahertz Oscillations in Graphene: Crater and Terraces in the Carrier Distribution Function

In graphene, after the electric field is turned-on, the ballistic acceleration of charge carriers up to the monochromatic optic phonon energy generates a back-and-forth motion of the whole distribution function between the zero point energy and the phonon energy. This effect is predicted to manifest in damped terahertz oscillations of the carrier drift velocity and average energy. The most dramatic feature of this transient phenomenon is the onset of momentum-free areas surrounded by high momentum probability in phase space, and smooth steps or terraces in the distribution function. This dynamical effect that only takes place within a voltage and sample length window, is the direct consequence of the interplay between the electric force and the randomizing nature of deformation potential optic phonons in the linear band structure of graphene.

Since its isolation in monolayer sheets, 1 graphene has emerged as a new and ideal two-dimensional (2D) material with exotic electronic properties. 2 The particular nature of its band structure with a linear dispersion between energy and momentum E =hv F |k|, where v F ∼ 10 8 cm/s is the Fermi velocity of electron and holes, results in a zero-gap material, where the overlap between the valence and conduction band is reduced to a single point around the K and K' symmetry points of the Brillouin zone. 2,3 Hence, charge carrier dynamics can be described by a formalism similar to Dirac relativistic equation, where the speed of light is replaced by v F which provides new opportunity to observe "quantum electro-dynamic" effects in nanoscale systems. For electronics applications, graphene is similarly attractive because of the linear energy-momentum dispersion, for which all charge carriers move with the velocity, v F much larger than in conventional semiconductor materials, 4 which anticipates faster time response.
Indeed, the carrier mean free path for collisions with low energy acoustic phonons (AP) can reach several micrometers, while efficient scattering with monochromatic optic phonons (OP) occurs at much larger energy (hω op ∼ 0.2eV ) 5 than in conventional semiconductors (hω op ∼ 0.04eV ). 6 Therefore, in electric fields F high enough (to escape AP or low energy scattering), charge carriers can experience quasi-ballistic run-away until they scattered with efficient OP emission once E ≥hω op . 7 The coherent aspect of this process i.e. the coherent acceleration of the carrier distribution function followed by quasi-instantenous carrier relaxation by OPs at high energy is expected to produce oscillations of the carrier velocity with a periodicity given by τ gr = hω op /eFv F . 8 In this paper we investigate carrier dynamics in graphene single layers, once the electric field has been turned-on, to determine the conditions of occurrence of current oscillations. Indeed, the thermal broadening of the initial distribution introduces a decoherence in the OP relaxation amongst carriers, which asides from low energy scattering, causes inherent damping of the oscillations, especially at high fields, where carriers overshoot the OP energy before emission. For this purpose, we use the Boltzmann formalism, and solve for the time-varying carrier distribution in the presence of OP scattering. [8][9][10] We also account for low energy scattering such as impurities and acoustic phonons, by using the relaxation time approximation. In weak concentrations (n c 10 11 cm −2 ), electron-electron interactions do not play a major role in transport in graphene, 11 and are not included in this analysis. We provide an analysis of the interplay between applied electric fields and strength of the low energy scattering rates for the onset of oscillations.
Model. Let us consider a system of electrons in the graphene conduction band under the influence of an electric field along the x-direction. We divide the momentum space into two regions, I (k < k c ) and II (k > k c ) separated by a circle of critical momentum k c = ω op /v F , that corresponds to the electron kinetic energy equal tohω op (Fig. 1a).
In region I electrons undergo quasi-ballistic acceleration and weak scattering by low energy mechanisms (e.g. impurities or APs) until they reach region II where they lose their energy by OP emission, and scatter back to region I. In our model, the electric field is assumed to be low enough so that electrons are scattered efficiently from region II to region I by OP emission with little probability to reach E ≥ 2hω op .
In the two regions, the time-dependent Boltzmann equations read: where f I ( k,t) and f II ( k,t) are the time-dependent momentum distribution functions in regions I and II respectively, and F is the electric field applied. ? The first terms in RHS of equations (1) account for low energy scattering mechanisms where, is the Fermi-Dirac equilibrium distribution function, and τ is the relaxation time (for the sake of simplicity, we assume τ is k-independent and we vary it's value compared to the OP scattering rate). S( k ′ , k) is the OP transition rate from a state with momentum k ′ to the state with momentum k, given by 12 where D o is the optical deformation potential, σ is the mass of the graphene sheet per unit area, and A is the area of the sheet. The second term on the right-hand side (RHS) of Eq. 1b is the carrier depopulation by OP emission, while in Eq. 1a it is the corresponding carrier repopulation at low energy. Here we neglect OP absorption processes sincehω op ≫ k B T , even at room temperature where the phonon occupation number N q is negligible, so OPs only scatter electrons from region II to region I.
We set the initial distribution (t = 0) as the Fermi-Dirac distribution, and solve Eqs.1 iteratively by noticing that for times 0 ≤ t <hk c /eF, the inside distribution drifts towards the critical circle with "speed" eF/h. During this trip, one assumes that all carriers are in region I, corresponding to f II ≈ 0. Eq. 1a is just a differential equation for which a solution is readily obtained. Then we solve Eq. 1b for k ≥ k c , which is of the same form as Eq. 1a but with the additional OP emission term. Here we use the boundary condition f I (k = k c ,t) = f II (k = k c ,t). 13 For t ≥hk c /eF, we substitute the solution f II ( k,t) into the integral of the right hand side of Eq. 1a to start the same procedure for later times.
As the electron population moves back and forth between regions I and II, each time undergoing more dephasing and broadening due to the finite duration of the OP emission process, we write the distribution function in each region at any time t corresponding to the n th trip toward the critical circle, i.e. for n = integer[eFt/hk c ] + 1, as a superposition of distributions f (i) and We note that the concentration n c in the conduction band is found to be n c ≈ 8 x 10 10 cm −2 for E F = 0 (at the Dirac point) used throughout this analysis. For Fermi levels at E F = k B T and E F = 2k B T above the Dirac point, the carrier concentrations would be n c ≈ 1.8 x 10 11 cm −2 , and n c ≈ 3.4 x 10 11 cm −2 , respectively. Results. We introduce the dimensionless time parameter β given by where t * =h k c eF is the approximate time taken for the distribution to complete one "trip". Fig. 2a shows the time evolution of the distribution function distribution (DF) for F = 1000V /cm (t * ≈ 2ps), in the absence of low energy scattering (γ = τ o /τ = 0). At β = 0(t = 0) the initial dis-tribution is the tail of the Fermi distribution in the conduction band (Eq. 12), which is Maxwellianlike, and drifts along the k x direction (opposite to the F-field) to reach the OP energy at β = 1 (second panel). At that time, a "hump" appears around k = 0 (not yet visible on the DF graph) as the front electrons reach region II, where they emit OP's and scatter back inside region I, as seen in the second panel of fig. 2b. As time progresses β = 1.5, the "hump" develops into a dimple, while the remaining high energy electrons from the first trip continue their drift in region II, where they experience strong OP depopulation. As seen in the third panel of fig. 2b, the corresponding repopulation rate at low energy exhibits a crater-like shape, which is due to the randomizing nature of the deformation potential OP scattering. Indeed, as schematically shown in fig.1b, after OP emission by high energy electrons, all | k − k c | -values become equiprobable, which forms a drifting circle of high repopulation rate, with increasing k-radius as first trip electrons penetrate deeper in region II. This crater-like feature of the repopulation rate is the primary cause of the dimple in the low energy distribution function, which are areas in k-space where electrons have low probability to scatter. As time progresses, the DF "dimple" evolves into a crater-like shape, ( fig.2.a, fourth and fifth panel), and the successive depopulation-repopulation OP processes overlapping at low energy with different amplitudes form also smooth terraces in the low energy tail of the distribution as it approaches steady state ( fig.2.a , fourth panel). These morphological effects in electron distribution are unique to graphene as a result its linear band structure and the interplay of the quasi-ballistic acceleration and relaxation by high energy monochromatic OP's.   3 shows snapshots of the distribution function at β = 2.5 for varying low energy relaxation time, expressed in terms of γ = τ 0 /τ. From the figure, one can see that the crater in the DF that occurs for γ = 0 progressively disappears as γ increases. Indeed low energy scattering in region I re-distributes charge carrier momenta around k = 0, especially in the crater center. For this reasons, the amplitude of the distribution in region I increases around k = 0. Also, the distribution amplitude decreases in region II (k > k c ) as carriers spend more time in region I (k ≤ k c ). Similarly the distribution recovers a streaming profile along the k x -direction as conventional semiconductors. 8,13 Nevertheless, even for strong damping (γ = 0.5), the distribution is characterized by a jagged profile, which contains the front and back ridges of the crater remnant still caused by the cumulative effects of the OP scattering for backward and forward carrier relaxation. The current density on the plane is given by, where φ is the angle between k and the k x axis. Fig. 4a shows the current density as a function of time (β ) for F = 1000V /cm and different low energy scattering rates. For weak scattering, the current density overshoots its steady state value through damped oscillations, as a result of the back and forth motion of the distribution function ( fig.1). The weaker the scattering, the higher the current overshoot. For strong low energy scattering, the current converges monotonically toward its steady state value without oscillations. 8 Notice the stronger the scattering, the lower the steady state current value. The insert shows the corresponding DF cross-sections at k x = 0.5 and β = 3 relative to the initial distribution at k x = 0. As low energy scattering increases, the DF becomes narrower, taking a "streaming" profile in the electric field direction. This is due to the fact that as the low energy scattering increases, fewer electrons reach the high energy ( E ≥hω op ) region II, reducing the number of electrons scattered back to the low energy region I, which narrows the distribution. Fig. 4b shows the current density as a function of time for three different field values, with γ = 0.05 in each case. One observes that the oscillation period τ gr =hω op /eFv F scales with the inverse of the electric field F. Also quite expectedly, the overshoot value increases with electric fields, but the damping is also enhanced with electric fields, which is due to the fact that the electron distribution penetrates the high energy(E ≥hω op ) region(II) farther thanhω op , stretching the carrier relaxation by OP emission, which in turn broadens the DF in region I along the field, thereby reaching steady state quicker. The insert shows the DF cross-sections at k x = 0.5 and β = 3 for the different fields, relative to the initial distribution at k x = 0. As the field is increased, the distribution function broadens, because more electrons reach the high-energy (E ≥hω op ) region II, and scatter back into region I broadening the distribution in the process.
The average energy density of carriers in the conduction band is given by, where f ( k,t) is the overall normalized distribution function, and the summation is taken over all k in regions I and II. At β = 0, the average energy per charge carrier is found to be E(0) ≈ 2.2k B T . This value is roughly twice larger than k B T expected for two-dimensional systems, and is the direct consequence of the linear energy dispersion in graphene, by contrast to the parabolic dispersion in normal 2D systems. Fig. 5a displays the ratio ε = E(β )/E(0) for F = 1kV /cm and different low energy scattering rate (γ). As expected, the energy converges to higher values as low energy scattering is reduced. In addition, the convergence toward steady state occurs through damped oscillations, even for γ = 0 as a consequence of the back and forth motion of the DF between the OP energy and the carrier zero-point energy. Quite clearly the oscillation period is given by t * for all γ. It is also seen that the oscillations persist even with significant low energy scattering. Fig. 5b displays the normalized energy ε as a function of time for different fields and γ = 0.1. As expected, carrier energies reach higher values as the field is increased, and the oscillation period decreases (larger β -period). Discussions. We have provided a transient analysis of the onset of current oscillations at the electric field turn-on caused by the back and forth motion of carrier distribution function between the zero-point energy and OPs in the presence of varying damping mechanisms. In this context we point out the anomalous shape of the carrier distribution as an interplay between ballistic acceleration and deformation potential OP emission in the transient regime. If OP-limited current os- cillations have been predicted in GaAs, 8 and indirectly observed in slightly n-doped InSb, 14 their manifestation in graphene is different in several respects: First, owing to the linear carrier energymomentum dispersion, the oscillation periodicity is given by τ gr =hω gr op /eFv F in graphene, while in GaAs parabolic conduction band it is expressed as τ GaAs = 2m * h ω GaAs op /eF, which nevertheless yields similar values, since the small effective mass, and the OP phonon frequency in GaAs compensate for the large v F . 8,10 Second, in III-V semiconductors, the OP polar nature focuses the low-energy repopulation along the k x axis, which provides a "streaming" profile to the carrier distribution instead of a crater-like shape in this case. Finally, in compound semiconductors, the oscillation onset is restricted by two conflicting conditions: On the one hand, the low value of the OP energy (hω op : 0.04eV ) is comparable to the thermal broadening of the carrier distribution at room temperature so that the back and forth motion of the distribution between the optic phonon and the zero point energy is immediately damped. [8][9][10] On the other hand, at low temperature, ionized impurity scattering becomes dominant and produces strong damping which can only be reduced by lowering the dopant density, thereby lowering the carrier density, and weakening the oscillation amplitude. In this respect, the high conductance of graphene, and the high optic phonon energy provide the conditions for room temperature observation.
Low energy scattering should however still be minimized. Usually, for experimental studies and device applications, graphene-layers rest on a dielectric substrate, or are confined between two dielectric slabs. 15,16 In this case the presence of high K-dielectrics sandwiching the strictly 2-D graphene-layer may be used to screen charged impurity that may reduce low energy elastic scattering. 17 However, dielectrics also contain interface and remote static charges that may offset dielectric screening. 18 Moreover, the interaction between the 2D carriers in graphene and low energy Remote Interface Phonon (RIP) arising from the proximity of the substrate 19,20 introduces new scattering sources. 21,22 Therefore in general, suspended graphene-layers avoiding the RIP influence may be preferable. 23 In this context to be observable at room temperature, the velocity oscillations should also take place within a parameter window. On the one hand the process requires eV >hω op where V is the external bias, to reach the OP energy, and on the other hand, the field should be large enough for carriers to escape low energy scattering, buthω op /eFv F τ op ≫ 1 where 1/τ op is the OP scattering rate, so that charge carriers don't penetrate the high-energy region E ≥hω op , scatter immediately after they reach the OP energy, which maintains the coherence of the distribution function. These requirements impose a lower and upper bound on the electric field, i.e. 0.5kV /cm < F ≪ 5kV /cm (in graphene), and a lower bound on the sample length L > v F τ op (> 1µm for τ op < 1ps), but L should be smaller than a few values of λ =hω op /eF, so as to prevent oscillation damping.
One important issue for the validity of this analysis is the effect of leakage current due to electrons in the valence band crossing to the conduction of carriers through the Dirac point. Indeed it has been shown that there is still a minimum conductance (G ∼ 4e 2 /h) between the two bands despite the singular nature of the Dirac point. 24 This value was later measured to be G ∼ e 2 /h. 1 However, this effect becomes important only when the graphene layer width w ≥ 23µm, for which the band gap E g ≤ 0.18meV are vanishing. Moreover, thermal effects such as acoustic phonon absorption by carriers in the valence band, are forbidden by conservation of both, energy and momentum. As for OP absorption we have shown earlier their occupation number is also negligible over the time scale considered in this analysis.

Supporting Information
Solution of Eqs. 1a & b. With the initial distribution (t = 0) given the Fermi-Dirac distribution, the solution for the "first trip" (0 ≤ t < k c /eF ) in region I (k ≤ k c ) is, where 1/ν o is the "speed" in k-space with which the distribution drifts in the field towards the critical circle.

High Energy Distribution Function
for k > k c even at room temperature, the approximation f o (k) ≈ 0 will be used for E ≥ ω op . By setting a new variable u = t − k x /eF , Boltzmann equation (1.b) in region II reads, and the total OP emmission rate in Eq. 1b (k > k c ) is, If we use ω op = 0.2eV , D o = 14eV /Å, σ = 7.61 x 10 −7 kg/m 2 , then τ o ≈ 0.32psec. The solution of Eq. S.3 is then, The solution (S.5) is valid only in the region −k c < k y < k c , which is the region of interest. Going back to the original variables, the general solution reads, The above solution describes a distribution function that drifts with "speed" eF/ in momentum space in the direction of the electric field, while decaying exponentially by OP emission.
Low Energy Distribution Function (k ≤ k c ). By using the same substitution as in the previous paragraph, the Boltzmann equation (1.a) inside the critical circle (region I) reads, and α(k) = arcsin{ k c k + k c } The solution to (S.7) is, in which the function Ω(k y ; u) is determined by initial conditions. Going back to the initial variables, By using the matching condition, The normalization condition determines χ(t) as Eqs. (S.11) do not exactly conserve the number of carriers n c . Indeed there is a decrease in the number of carriers in the system because at high fields electrons run away beyond region II i.e E ≥ 2 ω op without being scattered. Another reason is that the summation in Eqs. (S.11) is limited to the n = 3 terms thereby ignoring the contribution of later "trips"(n ≥ 4 terms) which might be relevant at higher fields considered in the analysis.