Nonperturbative model of harmonic generation in undoped graphene in the terahertz regime

We present a formalism to calculate the nonlinear terahertz response of monolayer graphene using a two-band tight-binding model of graphene. We develop the dynamic equations for the electron density matrix in the length gauge to calculate the interband and intraband carrier dynamics at terahertz frequencies. Using the calculated nonlinear interband and intraband current densities, we obtain the nonlinear transmitted and reflected terahertz fields. We find that when the conditions are such that the interband and intraband current densities are comparable, strong generation of odd harmonics is observed. We examine the response of undoped graphene as a function of a wide variety of system parameters, including operating temperature, Fermi velocity, pulse duration, and terahertz central frequency. We find that the response is generally most sensitive to the operating temperature and the terahertz central frequency. In particular, changing the ambient temperature from 10 to 100 K reduces the third harmonic field by a factor of 20, while increasing the terahertz central frequency from 0.5 to 5 THz results in an increase in the ratio of the generated third harmonic field to the fundamental of the reflected field from 19% to 49%. The very strong dependence on the temperature and central frequency is found to be due to the strong dependence of the nonlinear generation on the strength of the interplay between the intraband and interband dynamics. Our work demonstrates that, under the right conditions of low temperature and moderate THz field central frequency, third harmonic generation should be experimentally observable in undoped monolayer graphene at moderate THz field amplitudes of only 1 kV cm−1.


Introduction
Monolayer graphene has been shown to exhibit intriguing electrical, mechanical, and optical properties. Although it was first studied theoretically over a half century ago [1], it was brought to the limelight only in 2004 after its experimental discovery via mechanical exfoliation [2]. Many possible applications have been anticipated for graphene, such as photodetectors, transparent electrodes for displays, and terahertz (THz) transistors, to name just a few [3][4][5][6][7][8]. Applications in the THz domain are particularly interesting [9][10][11][12]. Many materials are highly transparent to terahertz radiation. This transparency, along with the relatively short THz wavelengths, offers an opportunity for high-resolution noninvasive imaging, sensing, high speed wireless communication and for security and quality-control applications [9][10][11]. However, there is a shortage of devices and components to efficiently generate, detect, and control terahertz waves. Recent efforts to generate single cycle intense THz pulses are opening a new avenue of exploration in condensed matter [13][14][15][16]. In particular, THz radiation is able to drive and probe electron transport within the bands in a wide variety of semiconductors and metals.
The zero-gap Dirac points along with the linear (as opposed to parabolic) electron band dispersion near the Dirac points [1,17,18], make the response of graphene to THz fields particularly interesting. The vanishing density of states where the two bands touch results in strong valence band to conduction band absorption for frequencies ranging from optical frequencies down to THz and below. The linear dispersion means that electrons and holes in graphene behave like massless 'relativistic' particles with an effective Fermi velocity » v c 300, F where c is the speed of light; this can result in a current density that depends nonlinearly on the driving THz field.
The absorption of optical and THz radiation by graphene can be characterized by interband and intraband transitions. The THz response is strongly dependent on the Fermi level, scattering mechanisms, ambient temperature, applied electric field amplitude, and central frequency. The low THz photon energy can prove very important in the study of intrinsic (undoped) graphene, where the chemical potential is at the Dirac point. For THz fields, interband transitions near the Dirac point take place, while the injected carriers are then driven within their respective bands. Earlier theoretical studies that employed the velocity gauge anticipated a strong nonlinearity in graphene at optical frequencies, and an even stronger one at lower frequencies such as terahertz [19][20][21][22][23].
Experimentally, third harmonic generation at terahertz frequencies has been observed using a 45-layer sample [24]. Yet no indication of third harmonic generation has been observed in monolayer graphene [25]. Almost all experiments have been performed using quite heavily doped samples, where the response is dominated by intraband dynamics. In this scenario, the clipping of the intraband current should lead to the generation of odd harmonics. However, if the carrier scattering is fast enough, the coherence of this process will be destroyed, which leads to diminished harmonics. For moderate Fermi energies that are a few tens of meV from the Dirac point, there can be a small contribution due to interband dynamics that is dependent on the field amplitude. If it is present, it leads to an increase in the current density, which reduces the current clipping and in turn reduces the harmonics. In our recent work [26] we showed that, if the Fermi energy is only within a few meV of the Dirac point, not only is the interband current comparable to the intraband one, but a strong nonlinearity in the interband current density can arise. In that work, we employed a two-band model in the length-gauge, considering suspended monolayer undoped graphene at a temperature of 10 K with an incident single-cycle THz pulse centered at a frequency of 1 THz. We showed that multiple harmonics of the THz field should be generated, and that the reflected third harmonic can be as large as 32% of the reflected fundamental. We did not, however, present a systematic study of the dependence of the nonlinearity on the system parameters.
In this paper, we present a detailed derivation of our two-band tight-binding model for the intraband and interband dynamics of the linear and nonlinear THz response of monolayer undoped suspended graphene. We then use this model to explore the sensitivity of the nonlinear response to a variety of key parameters. In particular, we investigate the role of the temperature, Fermi velocity, pulse duration, and central THz frequency on harmonic generation. The current densities and the corresponding harmonics are numerically calculated for two values of the temperature, 10 and 100 K. The results show a strong dependence of the third harmonic level on the operating temperature; for a central frequency of 1 THz, the ratio of the amplitude of the third harmonic to the fundamental in the reflected field drops from 32% at 10 K to about 1.6% at 100 K. More interestingly, we find that the ratio of the third harmonic to the fundamental initially increases with the incident field amplitude, but then peaks and decreases beyond a certain incident field amplitude. This arises due to the contribution of a fifth (and higher) order process to the generation of the third harmonic. Furthermore, an investigation of the magnitude of the Fermi velocity on the response shows that the peak value of the third harmonic is essentially independent of the Fermi velocity, but occurs at different field amplitudes for different Fermi velocities. In contrast, we find that, due to the effect of an increase in carrier injection, the peak of the third harmonic level decreases when the pulse duration is increased. Finally, and perhaps most importantly, we find that increasing the central frequency from 0.5 to 5 THz results in an increase by almost a factor of three in the ratio of the third harmonic to the fundamental in the reflected field.
The paper is organized as follows. In section 2, we present the derivation of our theoretical model. We start by using the tight-binding method to obtain interaction matrix elements in the length gauge, and then use these in a two-band density matrix formalism to derive the interband and intraband current densities and the transmitted and reflected fields. In section 3, we present the nonlinear THz response of numerical simulations for undoped graphene for different temperatures, Fermi velocities, pulse durations, and central THz frequencies. The conclusions are presented in section 4.

Theory
In order to calculate the nonlinear THz response of graphene, we first need to calculate the interband and intraband current densities. We accomplish this by employing a density matrix formalism in the basis of the graphene valence and conduction band states. Thus, before describing the density matrix approach, we first review the basic properties of the electron bands in graphene in the tight-binding approximation.

Tight-binding model
Graphene is formed from a single layer of carbon atoms arranged in a honeycomb lattice. We describe the hexagonal Bravais lattice using the primitive translation vectors a  a  a  a  x  y a  x  y  3  2   3  2  ,  3  2   3  2  .  1   o  o  o  o  1  2 The basis vectors, which form the A and B sublatices, are chosen to be = r 0, is the nearest-neighbor separation. The electronic states that are important for the THz response of graphene arise from the carbon 2p z orbitals, j ( ) r , pz which form the extended π-bands of graphene. In the nearest-neighbor tight-binding model for the π-electrons, the Bloch states formed from linear combinations of 2p z orbitals are given by: where W c is the area of a unit cell, the sum is over the Bravais lattice vectors, R and = { } n c v , labels the conduction and valence bands such as s = 1 n for the valence band and s = -1 n for the conduction band. The phase between the two sublattice states is given by k a k a i i 1 2 The states are normalized over an infinite volume so that The dispersion relations for the two bands are given by where g  3.033 eV is the nearest-neighbor coupling constant and E pz is the energy of the 2p states in carbon. As a result of the symmetry between the sublatices, the conduction and valence states in graphene are degenerate at the two so-called Dirac points: It is easily shown that for energies within a few hundred meV of the Dirac points, the dispersion is essentially linear. Thus, if we let then for small k, the band dispersion near both Dirac points becomes approximately where H o is the full Hamiltonian of unperturbed graphene, r is the electron position vector, = -| | e e is the charge of an electron and ( ) t E is the applied THz electric field. The carrier dynamics are calculated by solving the Heisenberg equation of motion for the density matrix in the basis of conduction band and valence band Bloch states. To solve these equations, we need the matrix elements of the Hamiltonian, which are The matrix elements of r between Bloch states can be shown to be given by [27,28]: where the dipole connection matrix elements are given by where ( ) u r nk is the periodic part of the Bloch function, which is given by We have evaluated these connection matrix elements using our tight-binding wavefunction. Ignoring the overlap of atomic wavefunctions on different atoms, we obtain and, when we restrict ourselves to the regions of k-space close to the Dirac point, we can use equation (16) for cos is the angular unit vector in cylindrical coordinates with the origin at the relevant Dirac point. Noting that these gradient terms diverge at the Dirac points, it is clear that they will dominate over the second term in the equation for x ( ) k .

nm
In fact, it is easy to show that to be consistent with our approximation for c d + ( ) K k , the final term in equation (23) should be ignored. Using this, we finally obtain for the intraband connection elements Similarly, for the interband connection elements, we obtain 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 disconnected, and calculate the dynamics around each of the Dirac points separately. In this approximation and staying within the linear dispersion regime, it can be shown that the current density coming from carriers around the two Dirac points are identical. Thus, in all that follows we shall restrict ourselves to carriers of one spin type that are only near the K Dirac point. Then in calculating and sums over carriers we shall include the factors of = g 2 which are the spin and valley degeneracies, respectively.
To simplify the following discussion and the calculation of the current densities, we will assume in all that follows that k is discretized on a uniform two-dimensional grid, such that Dirac delta functions in k become Kronecker delta functions, integrals become sums, and derivatives with respect to k are to be interpreted using finite difference approximations to the derivatives.
We define the matrix element of the density operator between two Bloch states to be Then, using the above expression for the matrix elements for the Hamiltonian, the dynamical equations for the density matrix elements become: is the Fermi-Dirac distribution with energy E, Fermi level m , F and a time-dependent temperature, T. We examine the dynamics of the vacancy population, r ( ) k , hh rather than the valence band population because the vacancies are only found close to the Dirac points, which enables us to just include states in relatively small regions in k-space about the Dirac points in our simulations. In the above equations, we have introduced phenomenological scattering terms to account for defect scattering, electron-phonon scattering and carrier-carrier scattering. For the interband coherences, r ( ) k , cv we introduce an interband decoherence time, t ( ) k and we force the populations to relax back to Fermi-Dirac quasithermal distributions, with relaxation times t ( ) k c and t ( ) k h for the electrons and holes, respectively. At each time step, a new temperature is chosen to obtain the carrier populations, which are time dependent due to the THz-induced interband transitions. Work by other researchers [29] has indicated that the time taken for conduction band electrons to relax to the valence band is much longer than intraband scattering times; we therefore neglect interband relaxation.
In all that follows, we will assume that at t = 0, the system is initially in thermal equilibrium before the THz pulse arrives. Thus, initial conditions for the diagonal density matrix elements are determined by the Fermi level, and the temperature of the system. At t = 0, before the THz pulse arises, there will be no off-diagonal elements in the density matrix. Hence, the initial condition will be: In equations (32) and (33), is the carrier quasithermal distributions that have a quasi-Fermi energy (and potentially temperature) that depends on time. These quasi-Fermi energies are determined at each time by requiring that where A is the the normalization area of the graphene sheet and s v hh k are the time-dependent electron and hole densities that one determines from the simulation at the given time. Note that the time dependence of the carrier densities (and hence temperature and quasi-Fermi energy) arises due to the THz-induced interband transitions.

Current density
In order to calculate the current density as a function of time, we follow the approach of Aversa and Sipe [27]. In order to be able to separate the calculations of the intraband and the interband current density, we split the position operator r into r i that has only diagonal components and r e that has only off-diagonal components in n.
Using this, the current density is given by where p is the electron momentum operator. The sums in the above equation are only over states in the vicinity of the K Dirac point and the factors g s and g v in all terms account for the Dirac-point degeneracy and the two spin states. 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 Using the density matrix dynamical equation, along with the cyclic property of the trace, we can write the above expression as is the interband contribution to the current and is the intraband contribution. The detailed steps for deriving the final equations of the current density are given in appendix A. Taking the exciting THz electric field to be x y the final expressions respectively for the total interband and intraband current densities are found to be: respectively. Note that the contributions given on the second and third lines for J i only contribute to third order in the electric field and higher. We generally find that they are small, but can be non-zero once the field introduces an asymmetry in r d + ( ) K k cv due to motion of the carriers in k-space.

Calculating the transmitted and the reflected fields
We take the graphene to be on a substrate that has an index of refraction n, assumed independent of frequency and that there is air just above the graphene. The THz field is incident normally on the graphene from the air and at the plane of the graphene is given by i i x i y Then, satisfying the boundary conditions (with our surface current), proceeding in the usual way [30] we find that the transmitted and the reflected fields at the position of the graphene are given by is the total graphene current density calculated using the transmitted field as the driving field and Z 0 is the impedance of free space. The details of how these equations are implemented self-consistently are given in appendix B.
There are a number of different ways in which one might try to solve the above dynamic equations. In this work, 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 introduce balanced difference quotient approximations [31] 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. It is quite easy to solve our dynamic equations for the density matrix components analytically to first order in the THz field. Doing this, and using our expressions for the interband and intraband current densities, leads to the standard equations for the interband and intraband conductivities. Furthermore, we have solved the full equations numerically and have verified that we obtain these expected linear conductivities. Finally, we have confirmed convergence in the nonlinear regime by changing the grid shape, grid density, the extent of the grid, the time-step tolerance and the polarization of the incident field.

Nonlinear THz response
In this section, we present the results of our simulations of the nonlinear response of undoped suspended graphene to THz fields. The input pulse is a single-cycle sinusoidal pulse with a Gaussian envelope as shown in figure 1(a), and is given by the following equation: where A o is chosen so that the maximum of the field, E o , is at the desired value 4 , t o = 3 ps is the time offset of the pulse, T FWHM is the full width at half maximum of the Gaussian envelope function and f o is the central frequency.
We have chosen to model undoped samples for two reasons. First, because most experiments and simulations to date have been performed on doped samples, there is considerably less understanding of the THz nonlinearity in undoped graphene. Second, we have found [26] that the ratio of the third harmonic to the fundamental in the current density generated in undoped graphene operating at low temperature is much larger than for doped graphene at the same incident field amplitudes. Therefore, the undoped system is very promising system for third harmonic generation.
In our previous work [26], we examined the effects of the scattering time on the nonlinear response. Different experimental studies have found a wide range of values for the scattering time that depend on the temperature, the doping level, the field amplitude and the substrate [32][33][34][35][36][37][38]. In this work, we only consider the response for fixed scattering and decay time constants of t t t = = = 50 fs c h . The choice of 50 fs is a conservative value for undoped graphene at low temperature [24,25,29,[39][40][41]. As we have shown [26], increasing this time will lead to stronger harmonic generation. Unless otherwise stated, in all of the simulations, we use the following parameters: temperature T = 10 K, Fermi velocityṽ c 300 F , THz pulse duration

Nonlinear response and the effect of temperature
We first present the current densities and reflected field spectra for different incident field amplitudes, and show the effect of raising the temperature from 10 to 100 K. As we shall see, the low temperature simulations give a much higher nonlinear response due to the enhanced interplay between inter-and intraband dynamics (schematically shown in figure 1(b)). In what follows, we examine both current densities and the reflection and transmission spectra, as shown in figures 2 and 3 respectively. Figure 2 shows the intraband (first row) and interband (second row) current density at 10 K (first column) and 100 K (second column) at four different THz field amplitudes, E o , of 1, 50, 100 and 200 V cm −1 . All the currents have been normalized to the corresponding peak amplitude, E o , of the incident field in order to facilitate a direct comparison. At 10 K, the normalized intraband current in figure 2(a) increases by a factor of three when the amplitude of the incident field is raised from 1 to 200 V cm −1 . This is expected as the low initial carrier density makes it possible for the applied THz field to inject carriers into the two bands, thereby increasing the intraband current density. The increase in the carrier density results in the filling of some originally open states in k-space, which introduces increased, time-dependent, Pauli blocking of the interband transitions. This in turn affects the interband current as shown in figure 2(b). As a result, the normalized interband current decreases and drastically changes form as the field amplitude is increased. This is at the heart of the predicted nonlinearity. We note that this mechanism takes place in a sub-cycle time scale and is beyond the conventional perturbative model, and full non-perturbative calculations such as we perform are necessary to declare the dynamics.
The transmitted and reflected spectral responses for T = 10 K are shown in figures 3(a) and (b), respectively. While there is no harmonic signal for the low field amplitude of 1 V cm −1 , the odd harmonics emerge as the field increases up to 100 V cm −1 . However, beyond that, the third harmonic begins to decrease. The third harmonic peaks at 0.06% of the fundamental of the transmitted signal, and remarkably at 32% of the fundamental of the reflected signal. This is due to the fact that the transmitted signal is dominated by the first term, is zero for suspended graphene, where n = 1. It is for this reason that we will now focus on reflected signal for suspended graphene. If a transmission geometry is still preferable for experimental reasons, a differential time-domain spectroscopy technique can be employed instead of the reflection geometry in order to measure essentially the same response [42,43].
We now consider the response for the same system when the temperature is raised. At T = 100 K, an order of magnitude increase in the intraband current is observed compared to the current at T = 10 K at a low field amplitude of 1 V cm −1 as shown in figure 2(c). This is simply a consequence of the larger number of thermal carriers in the two bands at the higher temperature. We also note that the increase in the relative intraband current with increasing field is essentially negligible at this temperature due to Pauli-blocking; the larger thermal population of carriers in this case limits the interband current, as the THz energy at 1 THz is just 4.14 meV, while = k T 8.6 B meV. This effect is also clearly seen in the interband current in figure 2(d), which is smaller than was found at T = 10 K. We also see that the change in the interband current with increasing THz field is much smaller than was found at T = 10 K. As a result of the Pauli-blocking, the interband contribution to harmonic generation is diminished. This is seen in figures 3(c) and (d), where for T = 100 K, the third harmonic amplitude is only 0.02% and 1.6% of the fundamental transmitted and reflected signal, respectively. The comparison of the response at the two temperatures emphasizes the role that the interplay between intraband and interband processes that is in the heart of the observed strong nonlinearity. It also clearly demonstrates the need to perform these experiments at low temperatures.

Effect of the Fermi velocity
We now turn to an examination of the effect of the Fermi velocity on the nonlinear response. In what we have presented thus far, we have taken the Fermi velocity to be the commonly chosen value of @´v 1 10 m s . However, a recent study [35] has shown that the Fermi velocity can range from0.85 10 6 to´-2.5 10 m s 6 1 , depending on the substrate. For simplicity we are considering suspended graphene, but it is generally important to understand if the magnitude of the Fermi velocity is critical to obtaining a strong harmonic response of graphene.
The Fermi velocity enters the simulations in two ways. First, it appears as a prefactor on the dominant term in the expression for the intraband current density (see equation (42)). Thus, an increase in the Fermi velocity results in an increase in the intraband current. However, because the main source of the third harmonic is the interband current, this increase in the intraband current does not significantly directly affect the third harmonic; however, it will affect the ratio of the third harmonic to the fundamental, due to the rise in the intraband component of the fundamental. The Fermi velocity also enters the interband dynamic equation through the energy difference, w ( ) k cv (equation (31)). Thus, an increase in the Fermi velocity will result in a shifting of carrier injection closer in k-space to the Dirac point, where the interband connection elements are larger. This, in turn, will result in a larger interband current density for a given electric field amplitude. For these reasons, we expect that an increase in the Fermi velocity will result in a more rapid increase in the third harmonic with electric field amplitude.
We present the results of simulations for three different values of the Fermi velocity in figure 4(a). Note that in all three cases the relative third harmonic amplitude initially increases, reaches a maximum, and then decreases. This general trend was noted in our previous work [26] and arises due to nonlinearities beyond third order. The dominant source of the third harmonic at low field amplitudes arises from the third-order nonlinearity, but as the incident field amplitude increases fifth order and higher responses produce a sub-cubic increase in the generated third harmonic field. The result is a clear maximum in the relative third harmonic response for all three different Fermi velocities. However, as expected from the discussion above, the required electric field to achieve the maximum third harmonic level decreases as the Fermi velocity is increased. This demonstrates that the results we have presented in previous figures for the maximum ratio of the third harmonic amplitude to the fundamental are robust with respect to changes in the actual Fermi velocity. But, because the maximum in the relative third harmonic amplitude occurs at higher incident field amplitudes when the Fermi velocity is lower, the absolute amplitude of the maximum third harmonic is higher when the Fermi velocity is smaller. In order to demonstrate this, we apply a bandpass filter to the reflected signal and retrieve the third harmonic time-domain signal 5 . Such a signal is shown in the inset of figure 4(b) for the case of = v v .

Fo F
The multi-cycle signal is a clear signature of the third harmonic. Now we plot the peak field amplitude of this third harmonic signal as a function of the incident field for all three values of the Fermi velocity as shown in figure 4(b). The field amplitude increases monotonically in all cases. At high fields, it is higher for systems with smaller Fermi velocities. The main reason for this somewhat surprising result is that the higher-order nonlinearities are weaker when the Fermi velocity is smaller; as a result the third harmonic is less suppressed by higher-order effects. Thus, although the third harmonic rises more slowly when the Fermi velocity is low, it is eventually larger if one goes to high enough incident field amplitudes. 5 A band-pass filter has been applied to the spectrum that is centered at the third harmonic peak and extended up to the minima from both sides. Then, an inverse Fourier transform has been applied in order to get the temporal third harmonic signal.

Effect of pulse duration and central frequency
Next, we study the effect of the THz pulse duration on third harmonic generation. In figure 5(a) we plot the third harmonic amplitude as a function of the THz field amplitude for pulse durations of T FWHM = 1, 1.5 and 2 ps. While there is almost no effect of the pulse duration on the response at low THz fields up to 50 V cm −1 , when the field is increased further we observe a decrease in the third harmonic ratio with increasing pulse duration. This arises because a longer pulse injects a larger number of carriers. This increases the carrier density in both bands.
Consequently, the intraband current is increased and, due to increased phase space filling, the amplitude of the interband current density (not shown here) is reduced. Because the interband current density is the main source of the nonlinearity, the reduction in the interband current reduces the nonlinearity. A similar trend is also observed in the absolute third harmonic amplitude, as shown in figure 5(b). Finally, we examine the effect of increasing the central frequency of the incident signal. The interplay between intraband and interband dynamics depends strongly on the operating central frequency, as can be expected from the simple picture shown in figure 6, where we see that a larger central frequency results in carriers being injected farther from the Dirac point. In figure 7(a), we plot the normalized third harmonic amplitude as a function of the incident THz field amplitude for central frequencies of 0.5, 2.0 and 5.0 THz. Because many strong THz sources are single-cycle pulses [10,14,44,45], we have adjusted the pulse durations in each case so that the product of the central frequency and the pulse duration is constant and the pulse remains single-cycle. The general trend of the third harmonic amplitude is similar to what has been presented above. However, compared to the above simulations at 1 THz, a larger maximum third harmonic amplitude is observed as the central frequency of the incident pulse is increased. For a THz central frequency of 2 THz, the normalized third harmonic amplitude peaks at about 46% at the relatively low incident field amplitude of 300 V cm −1 compared to a value of about 19% for the central frequency of 0.5 THz at a field amplitude of  75 V cm −1 . When the central frequency is increased to 5 THz, the maximum in the normalized third harmonic increases slightly to 49%, but the maximum is not achieved until the input field amplitude is 1 kV cm −1 .
Because the peak in the relative third harmonic signal occurs at higher incident field amplitudes for the higher-frequency pulses, the increase in the maximum of the amplitude of the generated third harmonic has an even stronger dependence on the central frequency than the normalized amplitude, as shown in figure 7(b). For a central frequency of 5 THz and incident peak amplitude of 1.2 kV cm −1 , the third harmonic amplitude is about 2.24 V cm −1 . We note that the absolute value of the reflected field at the fundamental decreases as we increase the central frequency, which leads to a difference in the values of the incident field at which the curves cross in figure 7(b) compared crossing points of the curves in figure 7(a).
We have performed simulations for frequencies up to 50 THz and find that the trends seen in figure 7 continue. First, the maximum in the normalized third harmonic amplitude continues to rise slightly with central frequency, but does not go much above 50%. Second, the higher frequencies have a lower nonlinearity at low field amplitudes but the nonlinearity peaks at a higher value when the central frequency is higher. This suggests that both the third order and fifth order response decrease as the central frequency of the pulse is increased. Therefore, if one wishes to obtain a strong third harmonic response, there is a trade-off. For a given input field, there is an optimum central frequency; for example for an input field of 100 V cm −1 , the optimum central frequency is approximately 1 THz. However, if the maximum field amplitude of the incident pulse is not an issue, then it appears that the response will be higher at higher central frequencies (up to at least 50 THz) for single-cycle pulses. The observed increase in the nonlinear interplay with increasing central frequency arises due to two factors. First, according to the Drude model for the linear intraband current density, we expect that the intraband current density will decrease with central frequency. Indeed, we find that at 200 V cm −1 the intraband current is suppressed by 19% as the central frequency is increased from 0.5 to 2 THz, as shown in figure 8(a). This decrease in the intraband current density increases the ratio of the third harmonic to the fundamental, simply due to the decrease in the reflected fundamental. Second, and more importantly, there is a significant increase in the interband current as the central frequency is increased. For instance, the interband current at 200 V cm −1 with a central frequency of 2 THz is more than double its value at a central frequency of 0.5 THz as shown in figure 8(b). This increase is largely arises because fewer of the interband transitions that are resonant with the 2 THz photons are occupied, either thermally or by injected carriers, as shown schematically in figure 6. Because almost all of the nonlinear response appears in the interband current, the larger interband current density results in a stronger interplay and hence a stronger third harmonic. We finally note that, as the central frequency increases, the THz pulse is shorter due to the fact that we have kept the ratio of the central frequency to the pulse duration constant. This leads to the injection of fewer carriers and consequently a significant increase in the field at which the third harmonic peaks as the THz frequency is increased, as is shown in figure 9.

Conclusion
In this work, we have presented a detailed study of harmonic generation in undoped suspended monolayer graphene at THz frequencies. A robust theoretical model based on the evolution of density matrix elements as determined by a length gauge technique has been developed. Expressions for the interband and intraband current densities and the transmitted and reflected THz fields were developed. We used simulations of the solutions of the dynamic equations to study the effects of various system parameters on the generation of the third harmonic, including temperature, Fermi velocity, pulse duration, central frequency of the incident pulse, and THz field amplitude. Our results show that the best conditions for harmonic generation are low temperature (∼10 K) and moderately high central frequency (∼5 THz). Under these conditions, we find that the third harmonic will be as much as 49% of the fundamental in the reflected signal, for a moderate incident THz field of only 1 kV cm −1 . Moreover, the absolute peak field increases monotonically with the amplitude of the incident field and reaches 2.24 V cm −1 (54.6 dB less than the transmitted fundamental) at incident field of 1.2 kV cm −1 . From practical point of view, either a transmission or reflection geometry can be employed to measure such a signal. However, if transmission geometry is preferable and the noise level is found to be a limiting factor, a differential time-domain spectroscopy technique can be employed instead in order to measure essentially the same response [42,43]. Moreover, a high sensitivity detection technique that allows a dynamic range larger than 60 dB is crucial in order to measure a 54.6 dB signal [46]. Such a value is indeed achievable, as it     From equations (26) and (27), we see that nn nn furthermore, it is easily shown that at all times, we must have r d

cc cc
It thus follows that when we sum over the entire first Brillouin zone, the expressions on the first two lines in equation (A.12) will sum to zero. Thus, we finally obtain If we restrict ourselves to regions in k-space near the two Dirac points, it is easily seen that the current density arising from sums over states near the two points are equal. Thus, we need only to perform calculations for the current density near the K-point. Using our expressions for the energy and connections elements near the K -point and multiplying by two for spin and additional two to account for the ¢ K -point gives us the final expression found in equation (42).

Appendix B. Field calculations
In this appendix, we present the details of the calculation of the transmitted field that must be employed in our numerical implementation of the density matrix equations. The key issue is that at each time step the evaluation of the derivative of the density matrix requires that we know the THz field at the graphene, which is the transmitted field. However, the calculation of the transmitted field at a given time requires that we know the current density at that time, which itself depends upon the transmitted field, as is seen in equation (44). We will show in the following that we can solve this problem by splitting the terms in the current equation into two parts that are dependent and independent of the transmitted field. We start by writing the current density in the form where « S is a matrix that is given by   Using this we finally obtain the following expressions for the transmitted and reflected fields at time t: