Abstract
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.