{'fraction_non_alphanumeric': 0.0359204618345093, 'fraction_numerical': 0.013761735378156162, 'mean_word_length': 4.566374553716326, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High harmonic generation in undoped graphene: Interplay of inter- and intraband dynamics

We develop a density matrix formalism in the length gauge to calculate the nonlinear response of intrinsic monolayer graphene at terahertz frequencies. Employing a tight-binding model, we find that the interplay of the interband and intraband dynamics leads to strong harmonic generation at moderate field amplitudes. In particular, we find that at low temperature, the reflected field of undoped suspended graphene exhibits a third harmonic amplitude that is 32% of the fundamental for an incident field of 100 V/cm. Moreover, we find that up to the seventh harmonic and beyond are generated.

As the paradigmatic example of a system exhibiting zero-gap Dirac points and linear band dispersion [1,2], graphene has been the subject of a host of studies, many with a focus on its fundamental physical and chemical properties, and many with a view towards device applications [3][4][5][6][7][8][9][10][11][12][13][14]. The absorption of electromagnetic radiation in graphene is controlled by interband and intraband transitions [2], as schematically shown in Fig. 1(a). In updoped graphene, applied fields lead to interband dynamics by inducing transitions between the bands, while the subsequent driving of the carriers within their bands by those same fields leads to intraband dynamics. Transitions near the Dirac point are accessed by terahertz (THz) fields, while for the intraband dynamics to result in significant currents the fields must have a relatively low frequency, as do THz waves [15][16][17].
Thus while the recent development of sources of intense THz pulses has led to new avenues of exploration across all of condensed matter physics [18][19][20][21], terahertz studies should be singularly suitable for revealing the physics of carrier dynamics in graphene. A typical THz single cycle pulse and its spectrum are shown in Figs. 1(b) and 1(c), respectively. Earlier theoretical studies anticipated strong optical nonlinearity in graphene, leading to effects such as high harmonic generation [22][23][24][25]. Although a recent experimental study has demonstrated generation of a third harmonic on the order of 10 −3 of the transmitted field power for a 45-layer graphene sample [26], there has been no experimental indication of such nonlinearities at terahertz frequencies using monolayer graphene [27].
In this paper, we develop a theoretical approach to treat the nonlinear interaction of THz radiation with monolayer graphene based on a density-matrix formalism in the length gauge. We predict that a strong third harmonic field will be emitted at incident fields as low as 100 V/cm, along with emission up to the 7 th harmonic and beyond, due to the interplay between the interband and intraband motion of the carriers. The critical requirements for such strong nonlinear effects are that the chemical potential is near the Dirac point and that the experiments are performed well below room temperature.
We start with a presentation of our theoretical model. We use a nearest-neighbor tight-binding model to treat the π−electrons in the graphene. The tight binding expression for the Bloch states is given by: where N is a normalization factor, n = {c, v} labels the conduction and valence bands and the sum is over the Bravais lattice vectors R. The factor, C = σ n e −iχ(k) , where σ v = 1 and σ c = −1 and χ (k) = arg [F (k)] is the phase between the two sublattice states, where F (k) ≡ 1 + e −ik·a1 + e −ik·a2 , with a 1 and a 2 being the primitive Bravais lattice translation vectors. The ϕ pz (r) are the 2p z orbitals of carbon, and δr A,B ≡ r − R − r A,B , where r A , r B are the basis vectors. In what follows, to simplify the notation we take the origin in k-space to be at K-Dirac point (where F (k) = 0). Near the Dirac point, it is easy to show that χ (k) θ + π/2 where θ is the angle that k makes with the k x -axis. In addition, the dispersion for the two bands is given approximately by E n (k) E pz − σ n v F k, where v F is Fermi velocity and E pz is the energy of the 2p z states in carbon.
It has been shown that the nonlinear response of semiconductors can only be reliably treated in a two band  [28,29]. In particular, if the more common velocity gauge is used, unphysical divergences arise in the nonlinear response at low frequencies, such as in the THz range, that can only be removed by developing sum rules; these become extremely complicated if one works to high order in the field [29]. Thus, in this work, we employ the length-gauge Hamiltonian, given by where H 0 is the full Hamiltonian of unperturbed graphene, e = −|e| is the charge of an electron, r is the electron position vector and E(t) is the THz electric field at the graphene. The carrier dynamics are calculated by solving the equations of motion for the reduced density matrix in the basis of conduction band and valence band Bloch states. These equations require the matrix elements of the Hamiltonian, which are given by The matrix elements of r between Bloch states can be shown to be given by [28,30]: where the connection elements, ξ nm (k), are given by where Ω c is the area of a unit cell and u nk (r) is the periodic part of the Bloch function. We have evaluated these connection elements using our tight-binding wavefunction. Ignoring the overlap of atomic wavefunctions on different atoms, near the Dirac point the connection elements are given approximately by Due to the large energy barriers between the two Dirac points and the small energies of the THz photons, we can treat the dynamics of the electrons near the two Dirac points as being disconnected and simply sum their (identical) contributions to the current density. We define the reduced density matrix elements to be ρ nm (k) ≡ a † m,k a n,k , where a † n,k (a n,k ) is the creation (annihilation) operator for an electron in the Bloch state |nk . Using Eq. (3) for the matrix elements of the Hamiltonian, the dynamic equations for the reduced density matrix elements become: and where ω cv (k) 2v F k and f n (k, t) is a Fermi-Dirac distribution with a time-dependent temperature. Note that the electric field in both equations is the field at the graphene, which is equal to the amplitude of the transmitted field. In our numerical implementation, we model the vacancy populations rather than the valence band electrons to allow us to only include states near the Dirac point. Defect scattering, electron-phonon scattering and carrier-carrier scattering are all potentially important processes in graphene. Because the scattering times in graphene are only on the order of tens of femtoseconds [26,27,[31][32][33][34][35], the inclusion of scattering processes is an essential element in any model of the THz response of graphene. In the above equations, we treat the scattering phenomenologically. For the interband coherences, ρ cv (k), we introduce an interband decoherence time, τ , which we assume is independent of k. The populations relax back to Fermi-Dirac thermal distributions, f n (k, t), with relaxation times, τ n and a temperature that is chosen to obtain the carrier populations, which are time-dependent due to the THz-induced interband transitions. As it has been found in previous studies [36] that the time taken for conduction band electrons to relax to the valence band is much longer than intraband scattering times, we neglect interband relaxation.
To solve the above equations, we employ a direct approach, where we put k on a grid and step through time using a Runge-Kutta algorithm. In order to facilitate this, we employ balanced difference approximations to the gradients. Given the geometry of the lattice and Brillouin zone, we employ a hexagonal grid with a uniform point density in k-space. Following the formalism of Aversa and Sipe [28], the current density is given by where the trace is over single-electron states, A is the normalization area of the graphene sheet, p is the electron momentum operator and ρ (t) is the reduced density matrix with matrix elements ρ nm (k). Using our expression for the Hamiltonian and the matrix elements of the position operator, it is possible to write the current density as the sum of an interband current density, J e , and an intraband current density, J i . After considerable work, one can show that the interband current density is given by [37] J e (t) = 8 |e| A Re while the intraband current density is given by Note that because to first order ρ cv (k) is symmetric under inversion, the terms on J i (t) given on the second line of Eq. (12) are to third order and higher in the electric field and are very small for the conditions considered. We consider a suspended graphene sample such as employed experimentally by Paul et al. [27] and use the time-dependent current sheet densities to calculate the transmitted and the reflected THz fields [38]. We now present the results of our simulations of the nonlinear THz response for undoped (µ F = 0) monolayer graphene. We run the simulation for a temperature of 10 K (corresponding to intrinsic density of 9.32 × 10 7 /cm 2 ) so as to minimize initial Pauli-blocking and thus maximize the interband transitions. The input THz pulse is a sinusoidal Gaussian pulse with central frequency of 1 THz and duration (FWHM) of 1 ps as shown in Fig.  1(b). The relaxation times τ and τ n have been chosen to be 50 fs, which is an average of various theoretical and measured values for such (sample-dependent) constants [26,27,[31][32][33][34][35]. This system has a thermal distribution of carriers with an average electron energy of 0.86 meV, which is much less than average photon energy of 4.14 meV of the 1 THz pulse. In the simulation, we keep the chemical potential at the Dirac point but as carriers are injected we raise the temperature in the distributions, f n (k, t), to which the carriers relax to account for the increase in the carrier density. All of the carriers (injected as well as thermal) are driven by the applied electric field (intraband dynamics). This results in the opening and closing of different states for extra injection of carriers. This interplay between intraband and interband dynamics is at the heart of the highly nonlinear response that we now examine.
In Fig. 2 we plot the intraband and interband current densities for four different incident field amplitudes (1, 50, 100 and 200 V/cm). In order to facilitate the comparison, all current densities are normalized to the peak, E i , of the incident field, such that, if the response were linear, these relative currents would be unchanged by an increase in incident field. To demonstrate the importance of the interplay of interband and intraband dynamics, we first present the results of the calculation when the interplay is missing. In the first case ( Fig. 2(a)), we have set the connection elements (ξ's) to zero in Eqs. (7) and (8) such that there is no electron motion between the bands. In this case, only the intraband current contributes. The relative current is decreased near the peaks (clipping) when the electric field is increased. This clipping results  7) and (8)  in a 23% reduction in the peak relative current when the field is 200 V/cm compared to the linear response; the reduction is expected due to the linear dispersion of graphene [1,2], but the nonlinearity is quite modest. In the second case ( Fig. 2(b)), all of the gradients in Eqs. (7) and (8) have been excluded. Interband transitions occur but the transfer of electrons between bands is not accompanied by the subsequent motion of electrons within individual bands. Here the relative interband current is seen to decrease as the THz electric field increases. This is due to interband absorption saturation that arises due to Pauli-blocking. Again, the nonlinearity is modest. Finally, we present the results of the full calculation that includes the full interband and intraband carrier dynamics. The relative intraband and interband current densities for this full calculation are shown in Figs. 2(c) and 2(d), respectively. The intraband current undergoes a large increase as the field is increased and is almost tripled when we go from a field of 1 V/cm to the highest field of 200 V/cm. This increase is expected, as it arises from the increase in the carrier densities due to the interband injection of carriers. In addition, there is a relatively small distortion in the intraband current as the field is increased. Finally, in Fig. 2(d), we plot the relative interband current density. We see first that as the incident field increases, there is a decrease in the relative interband current. Most importantly, the temporal form of the interband current is greatly modified, especially at fields at and above 100 V/cm. It is evident by comparison to Fig. 2(b) that this strong nonlinearity only arises when both interband an intraband processes are included in the calculation, and is thus due to the interplay between the interband and intraband dynamics. The normalized time-dependent reflected fields for different field amplitudes are shown for the full calculation in Fig. 3(a). The peak of the un-normalized reflected field is 0.34 V/cm for the incident field of 100 V/cm. As can be seen, the strong distortion of the interband current is clearly exhibited in the reflected field. This distortion is an indication of harmonics; the spectral responses normalized to the peak amplitude at the fundamental frequency are presented in Fig. 3(b). While there is no harmonic signal for the low field of 1 V/cm, the odd harmonics of the reflected field emerge as the field amplitude increases. The third harmonic peaks at 32% of the reflected spectral peak at the fundamental for a 100 V/cm incident field, and then decreases as the field is raised beyond that. Moreover, higher order odd harmonics such as the 5 th and the 7 th also appear as the field increases. We have also calculated the spectral response when there is no interplay between the interband and intraband dynamics and find (not shown) that the third harmonic amplitude is only 1.6% and 0.6%, respectively when the interband and intraband dynamics are shut off. Therefore, it is the interplay between the interband and intraband dynamics that is crucial to the appearance of the strong harmonic generation in undoped graphene. It is particularly noteworthy that we are obtaining a very strong nonlinearity even though the carrier densities are not very high at all; e.g., n = 2.24 × 10 8 /cm 2 for the 100 V/cm incident field. We have performed similar simulations under the same conditions but at room temperature, and find that the amplitude of the third harmonic generation is reduced by more than one order of magnitude. This reduction is due to the Pauli-blocking of the interband transitions due to the greatly increased thermal populations of carriers. We note that the amplitude of the generated third harmonic has a relatively weak dependence on the scattering times. For example, with a scattering time of 10 fs, the third harmonic for an incident field of 100 V/cm is still 10.5% of the fundamental.
Finally, to illustrate the combined effects of the interband and intraband transitions on the carrier distributions, we plot the electron density in k-space in Fig. 4 before the pulse arrives ( Fig. 4(a)) and at t = 2.75 ps, after almost half of the incident pulse has passed ( Fig.  4(b)). At t = 2.75 ps the distribution has changed in two key respects: (i) the carriers have been driven by the field so far to the right in k-space that most of them now have a velocity component in the positive-x direction, and (ii) the carrier density is considerably increased just above and below the k y = 0 line just to the right of the Dirac point as expected from Eq. (6). The first effect is the source of clipping in the intraband current density. The second effect results in the increase in carrier density, which resulted in the increased intraband current in Fig. 2(c). However, the most important effect of the strong redistribution of carriers in k-space is that it results in strong interplay between the interband and intraband dynamics. Although it is tempting to interpret this interplay as simply arising from the time-dependence of the Pauli-blocking of the states near the Dirac point, as a physical picture based on a rate equation model might suggest, such a picture does not lead to an accurate understanding of the detailed current dynamics because all of the dynamics are occurring on a sub-cycle timescale. Thus, the full simulation is required.
In conclusion, we have theoretically investigated nonlinear high harmonic generation in monolayer graphene. Our results demonstrate a very strong interplay between the intraband and interband dynamics, leading to large odd harmonics in the reflected field from suspended undoped graphene at low temperature. This work lays out the key conditions under which future experiments should be performed in order to achieve efficient high harmonic generation in monolayer graphene.
We thank the Natural Sciences and Engineering Re-