Third order optical nonlinearity of graphene

Due to the gapless linear dispersion of its low energy electronic excitations, graphene exhibits unique and remarkable optical properties: the linear optical response is characterized by a universal conductivity [1] ${{\sigma }_{0}}={{e}^{2}}/4\hbar$ at wavelengths from the mid-infrared to the visible, with a monolayer absorbance about $\sim 2.3%$ [2]. Thus while a single layer of graphene is highly transparent, considered as a monolayer it is a highly absorbing material, and saturated absorption is easy to achieve [3]. With its strong optical coupling, broadband absorption, and other novel material and electronic properties [4], graphene is a natural candidate for use in optically controlled devices in photonics and optoelectronics [5, 6]. In moving towards any application, an important step is understanding the optical nonlinearities of graphene. The graphene structure has center-of-inversion symmetry, and while second-order nonlinear response can arise from interface effects [7], nonuniformity of optical field [8], or the presence of dc currents [9], the first nonvanishing nonlinear susceptibility arising in pristine graphene is the third order susceptibility [10], and that is the focus of this paper.

Third order nonlinearities in graphene have been experimentally investigated; however, the extracted effective bulk susceptibilities $\chi _{\text{eff}}^{(3)}$ values show discrepancies and strongly depend on measurement method, light frequency, and perhaps sample preparation. In their pioneering work on the nonlinear response of graphene, Hendry et al [11] presented an expression for $\chi _{\text{eff}}^{(3)}$ for four-wave mixing experiments, and then approximated their expression to find a value of about ${{10}^{-15}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$ in the near-infrared, weakly dependent on frequency; they found agreement between experiment and their approximated expression by comparing their signal from graphene with that from a gold film. Using a sample with about 100 layers of graphene, Wu et al [12] found agreement with this approximated expression as well. However, in fact a direct calculation from the expression of Hendry et al gives a value only about ${{10}^{-19}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$ ³ . Yang et al [13] observed no detectable two-photon absorption in monolayer graphene by the Z-scan method at a wavelength of 780 nm, which limits $\text{Im}\left[ \chi _{\text{eff}}^{(3)} \right]$ to less than ${{10}^{-18}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$. The same experimental technique was used by Zhang et al [14] to measure the nonlinear refractive index at 1.55 μm, and they found ${{n}_{2}}\sim {{10}^{-11}}\;{{\text{m}}^{2}}\;{{\text{W}}^{-1}}$, six orders of magnitude larger than the value of $4.8\times {{10}^{-17}}\;{{\text{m}}^{2}}\;{{\text{W}}^{-1}}$ obtained at a similar wavelength by Gu et al [6], which matches a theoretical prediction using Hendry's original expression. Kumar et al [15] found a value of about $8\times {{10}^{-17}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$ by measuring third harmonic generation at a fundamental wavelength of 1720 nm, while Hong et al [16] found a value about two orders of magnitude smaller for third harmonic photon transitions at the M points in graphene. Sun et al [17] showed that two-color coherent injected currents can induce an observable terahertz radiation signal. Theoretically, besides the expression provided by Hendry et al, Rioux et al [18] investigated two-photon absorption and two-color coherent by using Fermi's Golden Rule for pristine graphene, and Jafari [19] calculated the nonlinear optical response in gapped graphene by adding a mass term. Zhang et al [20] used the density matrix method to study four wave mixing at the saturation regime in undoped graphene, and found an effective $\chi _{\text{eff}}^{(3)}$ about ${{10}^{-17}}$ m$^{2}$ ${{\text{V}}^{-2}}$ which decreased with increasing light intensity. A giant optical nonlinearity [21–23] in the presence of a strong magnetic field has also been calculated from a density matrix formalism. Third and higher harmonic generations from intraband contributions have been investigated in detail in the THz regime [10, 24, 25]. Yet a general calculation of even third order nonlinearities at higher frequencies, including the optical regime where interband contributions can play a central role, is still lacking.

A full calculation of the nonlinear response, even in the perturbative limit, would require the inclusion of electron–electron and electron–phonon scattering in the construction of a self-energy. Thermal effects caused by a high repetition rate of laser pulses, as used in Z-scan experiments, also cannot be ignored [26, 27]. In this paper, we neglect these complications, restrict ourselves to frequencies such that we can use approximations relevant around the Dirac cone, and perform a simple, zero-temperature perturbative calculation at the independent particle level, but fully taking into account both the interband and intraband motion of the electrons as perturbed by an incident field. The advantage of such a simplified approach is that we can obtain analytic results for the third order conductivity ${{\sigma }^{(3);dabc}}\left( {{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}} \right)$ in doped graphene; one can then determine an effective bulk susceptibility [11] by

$$\begin{eqnarray}&&\chi _{\text{eff}}^{(3);dabc}\left( {{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}} \right)={{\sigma }^{(3);dabc}}\left( {{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}} \right)/\left( -i{{\omega }_{t}}{{\epsilon }_{0}}{{d}_{\text{gr}}} \right),\end{eqnarray} \tag{ 1 }$$

with

$$\begin{eqnarray*}&&{{\omega }_{t}}={{\omega }_{1}}+{{\omega }_{2}}+{{\omega }_{3}},\end{eqnarray*}$$

${{\epsilon }_{0}}$ being the vacuum permittivity, and ${{d}_{\text{gr}}}\approx 3.3$ Å an effective thickness of graphene. These elementary but analytic results allow us to explore predictions for a number of third order nonlinear optical effects; they provide both an indication of what ranges of parameter space would be interesting to explore experimentally, and a benchmark for more sophisticated calculations.

For undoped graphene, we find that this model leads to a simple expression for the fully symmetrized conductivity

$$\begin{eqnarray}&&{{\sigma }^{(3);xxxx}}\left( {{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}} \right)=\frac{{{\sigma }_{0}}{{\left( \hbar {{v}_{F}}e \right)}^{2}}}{{{\hbar }^{4}}\left( {{\omega }_{1}}+{{\omega }_{2}} \right)\left( {{\omega }_{2}}+{{\omega }_{3}} \right)\left( {{\omega }_{3}}+{{\omega }_{1}} \right){{\omega }_{t}}}\end{eqnarray} \tag{ 2 }$$

where ${{v}_{F}}$ is the Fermi velocity. For doped graphene with a chemical potential $\mu$, however, our formulas show divergences related to the resonances between any involved photon energy and the chemical potential gap $2\left| \mu \right|$, with results quite different for different frequency combinations. Combined with the tunability of $\mu$ by an external gate voltage [28] or chemical doping [29] this should lead to novel approaches for controlling the nonlinear optical properties of graphene, and indeed to the possibility of graphene-based 'nonlinear optics on demand'.

We describe the electronic states in the $\pi$ and ${{\pi }^{*}}$ bands of graphene by a tight binding model employing carbon $2{{p}_{z}}$ orbitals; we denote by $\phi (\boldsymbol{r})$ centered at the origin. The Bloch states can be written as

$$\begin{eqnarray}&&{{\psi }_{s\boldsymbol{k}}}(\boldsymbol{r})=c_{A\boldsymbol{k}}^{s}{{\Phi }_{A\boldsymbol{k}}}(\boldsymbol{r})+c_{B\boldsymbol{k}}^{s}{{\Phi }_{B\boldsymbol{k}}}(\boldsymbol{r}).\end{eqnarray} \tag{ 3 }$$

Here $s\;=\;+\left( - \right)$ are band indexes for ${{\pi }^{*}}\left( \pi \right)$ bands, and ${{\Phi }_{\alpha \boldsymbol{k}}}(\boldsymbol{r})={{\mathcal{N}}^{-1}}{{\mathop{\sum }^{}}_{nm}}{{e}^{i\boldsymbol{k}\cdot {{\boldsymbol{r}}_{nm}}}}\phi \left( \boldsymbol{r}-{{\boldsymbol{r}}_{nm}}-{{\tau }_{\alpha }} \right)$ are tight binding basis functions for carbon atoms at sites $\alpha =A,B$ in unit cell ${{\boldsymbol{r}}_{nm}}=n{{\boldsymbol{a}}_{1}}+m{{\boldsymbol{a}}_{2}}$, where n and m are integers; $\mathcal{N}$ is the number of unit cells in the normalization area. We take the primitive lattice vectors to be ${{\boldsymbol{a}}_{1}}={{a}_{0}}\left( \frac{\sqrt{3}}{2}\hat{\boldsymbol{x}}\,-\frac{1}{2}\hat{\boldsymbol{y}}\, \right)$ and ${{\boldsymbol{a}}_{2}}={{a}_{0}}\left( \frac{\sqrt{3}}{2}\hat{\boldsymbol{x}}\,+\frac{1}{2}\hat{\boldsymbol{y}}\, \right),$ with the lattice constant ${{a}_{0}}=2.46$ Å. In our tight binding model, we set the energy reference by taking the onsite energy to be zero; we take the nearest neighbor coupling to be ${{\gamma }_{0}}=2.7$ eV [4], and denote the nearest neighbor overlap of the $2{{p}_{z}}$ orbitals as ${{d}_{0}}$. Neglecting overlaps of neighboring $2{{p}_{z}}$ beyond the first nearest neighbor, the Schrödinger equation for the $c_{\alpha \boldsymbol{k}}^{s}$ is

$$\begin{eqnarray}{{\gamma }_{0}}\left( \begin{array}{} 0 & {{f}_{\boldsymbol{k}}} \\ f_{\boldsymbol{k}}^{*} & 0 \\ \end{array} \right)\left( \begin{array}{} c_{A\boldsymbol{k}}^{s} \\ c_{B\boldsymbol{k}}^{s} \\ \end{array} \right)=\left( \begin{array}{} 1 & {{d}_{0}}\;{{f}_{\boldsymbol{k}}} \\ {{d}_{0}}\;f_{\boldsymbol{k}}^{*} & 1 \\ \end{array} \right)\left( \begin{array}{} c_{A\boldsymbol{k}}^{s} \\ c_{B\boldsymbol{k}}^{s} \\ \end{array} \right)\;.\end{eqnarray} \tag{ 4 }$$

Here ${{f}_{\boldsymbol{k}}}=1+{{e}^{-i\boldsymbol{k}\cdot {{\boldsymbol{a}}_{1}}}}+{{e}^{-i\boldsymbol{k}\cdot {{\boldsymbol{a}}_{2}}}}$ is the structure factor. The amplitudes of the eigenfunctions identified by equation (4) are then found to be

$$\begin{eqnarray*}&&\left( \begin{array}{} c_{A\boldsymbol{k}}^{s} \\ c_{B\boldsymbol{k}}^{s} \\ \end{array} \right)=\frac{1}{\sqrt{2(1+s{{d}_{0}}|{{f}_{\boldsymbol{k}}}|)}}\left( \begin{array}{} s\frac{{{f}_{\boldsymbol{k}}}}{|{{f}_{\boldsymbol{k}}}|} \\ 1 \\ \end{array} \right)\;,s=\pm ,\end{eqnarray*}$$

and eigenvalues are

$$\begin{eqnarray}&&{{\epsilon }_{s\boldsymbol{k}}}=\frac{s{{\gamma }_{0}}|{{f}_{\boldsymbol{k}}}|}{1+s{{d}_{0}}|{{f}_{\boldsymbol{k}}}|}.\end{eqnarray} \tag{ 5 }$$

Moving to a continuous range of crystal wavevectors $\boldsymbol{k}$, we write the Bloch functions as ${{\psi }_{s\boldsymbol{k}}}(\boldsymbol{r})=\frac{{{u}_{s\boldsymbol{k}}}(\boldsymbol{r})}{2\pi }{{e}^{i\boldsymbol{k}\cdot r}}$, where ${{u}_{s\boldsymbol{k}}}(\boldsymbol{r})={{u}_{s\boldsymbol{k}}}\left( \boldsymbol{r}+{{\boldsymbol{r}}_{nm}} \right)$ for any $n,m$. From the normalization of the Bloch functions $\mathop{\int }^{}d\boldsymbol{r}\ \psi _{s\boldsymbol{k}}^{*}(\boldsymbol{r}){{\psi }_{{{s}^{\prime }}{{\boldsymbol{k}}^{\prime }}}}(\boldsymbol{r})={{\delta }_{s{{s}^{\prime }}}}\delta \left( \boldsymbol{k}-{{\boldsymbol{k}}^{\prime }} \right)$ where the integral is over all space and we take $\boldsymbol{k}$ and ${{\boldsymbol{k}}^{\prime }}$ to be in the first Brillouin zone, the ${{u}_{s\boldsymbol{k}}}(\boldsymbol{r})$ are normalized according to

$$\begin{eqnarray*}&&{{\mathop{\int }^{}}_{\text{cell}}}\frac{d\boldsymbol{r}}{{{\mathcal{A}}_{\text{cell}}}}u_{s\boldsymbol{k}}^{*}(\boldsymbol{r}){{u}_{{{s}^{\prime }}\boldsymbol{k}}}(\boldsymbol{r})={{\delta }_{s{{s}^{\prime }}}},\end{eqnarray*}$$

where the indicated integral is over a unit cell with area ${{\mathcal{A}}_{\text{cell}}}$ in the xy plane and over all z. The Berry connections [30–32],

$$\begin{eqnarray}&&{{\xi }_{{{s}_{1}}{{s}_{2}}\boldsymbol{k}}}=i{{\mathop{\int }^{}}_{\text{cell}}}\frac{d\boldsymbol{r}}{{{\mathcal{A}}_{\text{cell}}}}u_{{{s}_{1}}\boldsymbol{k}}^{*}(\boldsymbol{r}){{\nabla }_{\boldsymbol{k}}}{{u}_{{{s}_{2}}\boldsymbol{k}}}(\boldsymbol{r}),\end{eqnarray} \tag{ 6 }$$

are then found to be

$$\begin{eqnarray}&&{{\xi }_{ss\boldsymbol{k}}}=\frac{1}{2}\left( \tau -{{\chi }_{\boldsymbol{k}}} \right),\;{{\xi }_{s\bar{s}\boldsymbol{k}}}=\frac{\tau +{{\chi }_{\boldsymbol{k}}}}{2\sqrt{1-d_{0}^{2}|{{f}_{\boldsymbol{k}}}{{|}^{2}}}},\end{eqnarray} \tag{ 7 }$$

with ${{\chi }_{\boldsymbol{k}}}=\text{Im}\left[ {{\nabla }_{\boldsymbol{k}}}{{f}_{\boldsymbol{k}}}/{{f}_{\boldsymbol{k}}} \right],$ $\tau ={{\tau }_{B}}-{{\tau }_{A}}$, and $\bar{s}$ the index of the band that is not the s band. At the Dirac points, ${{\chi }_{\boldsymbol{k}}}$ is singular.

Under the approximation that the optical field is treated as uniform, the light-matter interaction can in principle be described either through the use of a scalar potential or a vector potential associated with the electric field $\boldsymbol{E}\left( t \right)$. Both treatments have disadvantages. The first involves introducing the position operator $\boldsymbol{r}$, which does not have well-behaved matrix elements between the periodic Bloch functions ${{\psi }_{s\boldsymbol{k}}}(\boldsymbol{r})$. Starting with the work of Blount [30], a number of strategies and techniques have been developed for working with the effective matrix elements of the position operator that arise; these matrix elements are linked to the Berry connections [32]. Such problems do not arise if the vector potential is used to represent the electric field, but the inevitable band truncation in numerical calculations can result in false divergences [32, 33]. This problem does not plague calculations based on the position operator. And while these divergences can be eliminated using proper sum rules that involve all the bands, since only two bands are included in our tight binding model we use the approach based on introducing the position operator, and take the interaction Hamiltonian to be ${{H}_{I}}=-e\boldsymbol{E}\left( t \right)\cdot \boldsymbol{r}$, where we take the electronic charge to be $e=-\left| e \right|$. In the independent particle approximation we adopt in this paper, the proper treatment of the position operator then leads to the description of the system by the semiconductor Bloch equations [32]

$$\begin{eqnarray}&&\frac{\partial \rho _{\boldsymbol{k}}^{S}}{\partial t}=-i\left[ {{\mathcal{E}}_{\boldsymbol{k}}}-e\boldsymbol{E}\left( t \right)\cdot {{\xi }_{\boldsymbol{k}}},\rho _{\boldsymbol{k}}^{S} \right]-e\boldsymbol{E}\left( t \right)\cdot {{\nabla }_{\boldsymbol{k}}}\rho _{\boldsymbol{k}}^{S}.\end{eqnarray} \tag{ 8 }$$

Here $\rho _{\boldsymbol{k}}^{S}\left( t \right)$ is $2\times 2$ density matrix, with elements $\rho _{{{s}_{1}}{{s}_{2}}\boldsymbol{k}}^{S}\left( t \right)$, in which diagonal term $\rho _{ss\boldsymbol{k}}^{S}$ describes the population of band s at $\boldsymbol{k}$ and off-diagonal terms $\rho _{s\bar{s}\boldsymbol{k}}^{S}$ describe the optical polarization; similarly ${{\mathcal{E}}_{\boldsymbol{k}}}$ denotes the matrix with elements ${{\mathcal{E}}_{{{s}_{1}}{{s}_{2}}\boldsymbol{k}}}={{\delta }_{{{s}_{1}}{{s}_{2}}}}{{\epsilon }_{{{s}_{1}}\boldsymbol{k}}}$ and ${{\xi }_{\boldsymbol{k}}}$ denotes the matrix with elements ${{\xi }_{{{s}_{1}}{{s}_{2}}\boldsymbol{k}}}$. The areal current density is then calculated by

$$\begin{eqnarray*}&&\boldsymbol{J}=e\mathop{\int }^{}\frac{d\boldsymbol{k}}{{{(2\pi )}^{2}}}\text{Tr}[{{\boldsymbol{v}}_{\boldsymbol{k}}}\rho _{\boldsymbol{k}}^{s}],\end{eqnarray*}$$

where the integral is over the Brillouin zone and the ${{\boldsymbol{v}}_{\boldsymbol{k}}}$ are the matrix elements of the velocity operator $\boldsymbol{v}=-i{{\hbar }^{-1}}\left[ \boldsymbol{r},{{H}_{0}}+{{H}_{I}} \right],$ given by

$$\begin{eqnarray}&&{{\boldsymbol{v}}_{ss\boldsymbol{k}}}={{\nabla }_{\boldsymbol{k}}}{{\epsilon }_{s\boldsymbol{k}}},{{\boldsymbol{v}}_{s\bar{s}\boldsymbol{k}}}=i\left( {{\epsilon }_{s\boldsymbol{k}}}-{{\epsilon }_{\bar{s}\boldsymbol{k}}} \right){{\xi }_{s\bar{s}\boldsymbol{k}}}.\end{eqnarray} \tag{ 9 }$$

The direct perturbative solution of equation (8) has been discussed in detail before [32], in fact for a 3D crystal, using an approach in which interband and intraband contributions are separated as the perturbation expansion is developed. Here we introduce a simpler strategy by shifting to a moving frame that essentially follows the intraband motion of the carriers; we put $\rho _{\boldsymbol{k}}^{S}\left( t \right)={{\rho }_{\boldsymbol{k}+e\boldsymbol{a}\left( t \right)}}\left( t \right)$, where $\boldsymbol{a}\left( t \right)$ is the vector potential used to describe the electric field, $\boldsymbol{E}\left( t \right)=-\partial \boldsymbol{a}\left( t \right)/\partial t$. In a rough sense we are using the vector potential to treat the intraband motion, and the scalar potential to treat the interband motion. Equation (8) then reduces to

$$\begin{eqnarray}&&\frac{\partial {{\rho }_{\boldsymbol{k}}}\left( t \right)}{\partial t}=-i\left[ {{H}_{\boldsymbol{k}-e\boldsymbol{a}\left( t \right)}}\left( \boldsymbol{E}\left( t \right) \right),{{\rho }_{\boldsymbol{k}}}\left( t \right) \right],\end{eqnarray} \tag{ 10 }$$

with ${{H}_{\boldsymbol{k}}}\left( \boldsymbol{E}\left( t \right) \right)={{\mathcal{E}}_{\boldsymbol{k}}}-e\boldsymbol{E}\left( t \right)\cdot {{\xi }_{\boldsymbol{k}}}$. The ${{a}^{th}}$ Cartesian component of the areal current density is then written as ${{J}^{a}}\left( t \right)=e\mathop{\int }^{}\frac{d\boldsymbol{k}}{{{\left( 2\pi \right)}^{2}}}\text{Tr}\left[ v_{\boldsymbol{k}-e\boldsymbol{a}\left( t \right)}^{a}\;{{\rho }_{\boldsymbol{k}}}\left( t \right) \right]$.

We wish to generate a perturbative expansion of the solution of (10), which will then allow us to write an expansion of $\boldsymbol{J}\left( t \right)$ in terms of powers of the electric field,

$$\begin{eqnarray}&&\boldsymbol{J}\left( t \right)={{\boldsymbol{J}}^{(1)}}\left( t \right)+{{\boldsymbol{J}}^{(2)}}\left( t \right)+{{\boldsymbol{J}}^{(3)}}\left( t \right)+\cdots \;.\end{eqnarray} \tag{ 11 }$$

Because graphene has center-of-inversion symmetry the second order term ${{\boldsymbol{J}}^{(2)}}\left( t \right)$ will be identically zero in the dipole approximation.

We begin by writing

$$\begin{eqnarray}&&\boldsymbol{a}\left( t \right)=\mathop{\int }^{}\frac{d\omega }{2\pi }{{\boldsymbol{a}}_{\omega }}{{e}^{-i\omega t}},\end{eqnarray} \tag{ 12 }$$

and expand ${{H}_{\boldsymbol{k}-e\boldsymbol{a}\left( t \right)}}\left( \boldsymbol{E}\left( t \right) \right)$ in $\boldsymbol{E}\left( t \right)$ and $\boldsymbol{a}\left( t \right)$ to find

$$\begin{eqnarray}\begin{array}{} {{H}_{\boldsymbol{k}-e\boldsymbol{a}\left( t \right)}}\left( \boldsymbol{E}\left( t \right) \right) & = & {{\mathcal{E}}_{\boldsymbol{k}}}-\mathop{\int }^{}\frac{d\omega }{2\pi }eA_{\omega }^{a}{{e}^{-i\omega t}}\left[ \mathcal{V}_{\boldsymbol{k}}^{a}\left( \omega \right) \right. \\ {} & {} & \left. -\frac{e{{A}^{b}}\left( t \right)}{2}\mathcal{M}_{\boldsymbol{k}}^{ab}\left( \omega \right)+\frac{{{e}^{2}}{{A}^{b}}\left( t \right){{A}^{c}}\left( t \right)}{6}\mathcal{W}_{\boldsymbol{k}}^{abc}\left( \omega \right)+\cdots \right], \\ \end{array}\end{eqnarray} \tag{ 13 }$$

where superscripts indicate Cartesian components, to be summed over if repeated, and with $\mathcal{V}_{\boldsymbol{k}}^{a}\left( \omega \right)=\partial {{\mathcal{E}}_{\boldsymbol{k}}}/\partial {{k}^{a}}+i\omega \xi _{\boldsymbol{k}}^{a}$, $\mathcal{M}_{\boldsymbol{k}}^{ab}\left( \omega \right)=\partial \left[ \partial {{\mathcal{E}}_{\boldsymbol{k}}}/\partial {{k}^{a}}+2i\omega \xi _{\boldsymbol{k}}^{a} \right]/\partial {{k}^{b}}$, and $\mathcal{W}_{\boldsymbol{k}}^{abc}\left( \omega \right)={{\partial }^{2}}\left[ \partial {{\mathcal{E}}_{\boldsymbol{k}}}/\partial {{k}^{a}}+3i\omega \xi _{\boldsymbol{k}}^{a} \right]/\partial {{k}^{b}}\partial {{k}^{c}}$. In doped graphene, the singularities of ${{\boldsymbol{v}}_{\boldsymbol{k}}}$ and ${{\xi }_{\boldsymbol{k}}}$ at Dirac points can be ignored. For weak external electric fields $\boldsymbol{E}\left( t \right)$, equation (10) can then be solved iteratively. By expanding ${{\rho }_{\boldsymbol{k}}}=\rho _{\boldsymbol{k}}^{0}+{{\mathop{\sum }^{}}_{i}}\rho _{\boldsymbol{k}}^{\left( i \right)}\left( t \right)$ with $\rho _{\boldsymbol{k}}^{0}$ being the density matrix at equilibrium states and

$$\begin{eqnarray*}\begin{array}{} \rho _{\boldsymbol{k}}^{(1)}\left( t \right) & = & \mathop{\int }^{}\frac{d\omega }{2\pi }\left| e \right|A_{\omega }^{a}\mathcal{P}_{\boldsymbol{k}}^{(1);a}\left( \omega \right){{e}^{-i\omega t}}, \\ \rho _{\boldsymbol{k}}^{(2)}\left( t \right) & = & \mathop{\int }^{}\frac{d{{\omega }_{1}}d{{\omega }_{2}}}{{{\left( 2\pi \right)}^{2}}}{{\left| e \right|}^{2}}A_{{{\omega }_{1}}}^{a}A_{{{\omega }_{2}}}^{b}\mathcal{P}_{\boldsymbol{k}}^{(2);ab}\left( {{\omega }_{1}},{{\omega }_{2}} \right){{e}^{-i\left( {{\omega }_{1}}+{{\omega }_{2}} \right)t}}, \\ \rho _{\boldsymbol{k}}^{(3)}\left( t \right) & = & \mathop{\int }^{}\frac{d{{\omega }_{1}}d{{\omega }_{2}}d{{\omega }_{3}}}{{{\left( 2\pi \right)}^{3}}}{{\left| e \right|}^{3}}A_{{{\omega }_{1}}}^{a}A_{{{\omega }_{2}}}^{b}A_{{{\omega }_{3}}}^{c}\mathcal{P}_{\boldsymbol{k}}^{(3);abc}\left( {{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}} \right){{e}^{-i{{\omega }_{t}}t}}, \\ \end{array}\end{eqnarray*}$$

the iterative solutions are given as

$$\begin{eqnarray*}\begin{array}{} \mathcal{P}_{\boldsymbol{k}}^{(1);a} & = & {{I}^{\omega }}\left( \left[ \mathcal{V}_{\boldsymbol{k}}^{a},\rho _{\boldsymbol{k}}^{0} \right] \right), \\ \mathcal{P}_{\boldsymbol{k}}^{(2);ab} & = & {{I}^{{{\omega }_{1}}+{{\omega }_{2}}}}\left( \left[ \mathcal{V}_{\boldsymbol{k}}^{a},\mathcal{P}_{\boldsymbol{k}}^{(1);b} \right]+\frac{1}{2}\left[ \mathcal{M}_{\boldsymbol{k}}^{ab},\rho _{\boldsymbol{k}}^{0} \right] \right), \\ \mathcal{P}_{\boldsymbol{k}}^{(3);abc} & = & {{I}^{{{\omega }_{t}}}}\left( \left[ \mathcal{V}_{\boldsymbol{k}}^{a},\mathcal{P}_{\boldsymbol{k}}^{(2);bc} \right]+\frac{1}{2}\left[ \mathcal{M}_{\boldsymbol{k}}^{ab},\mathcal{P}_{\boldsymbol{k}}^{(1);c} \right]+\frac{1}{6}\left[ \mathcal{W}_{\boldsymbol{k}}^{abc},\rho _{\boldsymbol{k}}^{0} \right] \right). \\ \end{array}\end{eqnarray*}$$

Hereafter we keep the dependence of the terms on frequencies implicit, always linking the frequencies ${{\omega }_{1}}/{{\omega }_{2}}/{{\omega }_{3}}$ to the Cartesian components $a/b/c$. The matrix function ${{I}^{\omega }}\left( {{X}_{\boldsymbol{k}}} \right)$, where ${{X}_{\boldsymbol{k}}}$ is any $2\times 2$ matrix associated with wavevector $\boldsymbol{k}$, is defined as $I_{{{s}_{1}}{{s}_{2}}}^{\omega }\left( {{X}_{\boldsymbol{k}}} \right)={{X}_{{{s}_{1}}{{s}_{2}}\boldsymbol{k}}}/\left( \hbar \omega -\hbar {{\omega }_{{{s}_{1}}{{s}_{2}}\boldsymbol{k}}}+i\delta \right)$ with $\hbar {{\omega }_{{{s}_{1}}{{s}_{2}}\boldsymbol{k}}}={{\epsilon }_{{{s}_{1}}\boldsymbol{k}}}-{{\epsilon }_{{{s}_{2}}\boldsymbol{k}}}$. Expanding

$$\begin{eqnarray}&&v_{\boldsymbol{k}-e\boldsymbol{a}\left( t \right)}^{a}=v_{\boldsymbol{k}}^{a}-e{{A}^{b}}\left( t \right)M_{\boldsymbol{k}}^{ab}+\frac{1}{2}{{e}^{2}}{{A}^{b}}\left( t \right){{A}^{c}}\left( t \right)W_{\boldsymbol{k}}^{abc}-\frac{1}{6}{{e}^{3}}{{A}^{b}}\left( t \right){{A}^{c}}\left( t \right){{A}^{d}}\left( t \right)Z_{\boldsymbol{k}}^{abcd}+\cdots \;.\end{eqnarray} \tag{ 14 }$$

where $M_{\boldsymbol{k}}^{ab}=\partial v_{\boldsymbol{k}}^{a}/\partial {{k}^{b}}$, $W_{\boldsymbol{k}}^{abc}=\partial M_{\boldsymbol{k}}^{ab}/\partial {{k}^{c}}$, and $Z_{\boldsymbol{k}}^{abcd}=\partial W_{\boldsymbol{k}}^{abc}/\partial {{k}^{d}}$, we find that in an expansion of the ${{\boldsymbol{J}}^{(1)}}\left( t \right)$ of (11) according to (12) we have

$$\begin{eqnarray*}&&{{J}^{(1)d}}(t)=\mathop{\int }^{}\frac{d\omega }{2\pi }{{\sigma }^{(1);da}}(\omega ){{E}_{\omega }}{{e}^{-i\omega t}},\end{eqnarray*}$$

where the linear conductivity ${{\sigma }^{(1);da}}\left( \omega \right)$ is given by

$$\begin{eqnarray}\begin{array}{} {{\sigma }^{(1);da}}\left( \omega \right) & = & -\frac{{{e}^{2}}}{i\omega }\underset{\boldsymbol{k}}{\mathop{\mathop{\sum }^{}}}\,\text{Tr}\left[ v_{\boldsymbol{k}}^{d}\mathcal{P}_{\boldsymbol{k}}^{(1);a}+M_{\boldsymbol{k}}^{da}\rho _{\boldsymbol{k}}^{0} \right], \\ {} & = & \frac{i{{e}^{2}}}{{{\hbar }^{2}}\omega }\underset{s\boldsymbol{k}}{\mathop{\mathop{\sum }^{}}}\,M_{ss\boldsymbol{k}}^{da}{{n}_{s\boldsymbol{k}}}+\frac{{{e}^{2}}}{i\hbar \omega }\underset{s\boldsymbol{k}}{\mathop{\mathop{\sum }^{}}}\,\frac{v_{\bar{s}s\boldsymbol{k}}^{d}\mathcal{V}_{s\bar{s}\boldsymbol{k}}^{a}\left( \omega \right)s\left( {{n}_{+\boldsymbol{k}}}-{{n}_{-\boldsymbol{k}}} \right)}{\hbar \omega -s\hbar {{\omega }_{+-\boldsymbol{k}}}+i\delta }. \\ \end{array}\end{eqnarray} \tag{ 15 }$$

We also find

$$\begin{eqnarray}&&{{J}^{(3)d}}\left( t \right)=\mathop{\int }^{}\frac{d{{\omega }_{1}}d{{\omega }_{2}}d{{\omega }_{3}}}{{{\left( 2\pi \right)}^{3}}}{{\tilde{\sigma }\,}^{(3);dabc}}\left( {{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}} \right){{e}^{-i{{\omega }_{t}}t}},\end{eqnarray} \tag{ 16 }$$

where as in equation (1) we have ${{\omega }_{t}}={{\omega }_{1}}+{{\omega }_{2}}+{{\omega }_{3}}$, and the (nonsymmetrized) third order conductivity ${{\tilde{\sigma }\,}^{(3);dabc}}\left( {{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}} \right)$ is

$$\begin{eqnarray}&&{{\tilde{\sigma }\,}^{(3);dabc}}=\frac{-i{{e}^{4}}}{{{\omega }_{1}}{{\omega }_{2}}{{\omega }_{3}}}\underset{\boldsymbol{k}}{\mathop{\mathop{\sum }^{}}}\,\text{Tr}\left[ v_{\boldsymbol{k}}^{d}\mathcal{P}_{\boldsymbol{k}}^{(3);abc}+M_{\boldsymbol{k}}^{da}\mathcal{P}_{\boldsymbol{k}}^{(2);bc}+\frac{1}{2}W_{\boldsymbol{k}}^{dab}\mathcal{P}_{\boldsymbol{k}}^{(1);c}+\frac{1}{6}Z_{\boldsymbol{k}}^{dabc}\rho _{\boldsymbol{k}}^{0} \right].\end{eqnarray} \tag{ 17 }$$

In graphene, the hexagonal lattice has ${{D}_{6\,h}}$ symmetry, which results in zero second order conductivity. The linear conductivity has only one independent nonzero component ${{\sigma }^{(1);xx}}={{\sigma }^{(1);yy}}$. For the third order conductivity there are in all eight nonzero components, among which three are independent; we have

$$\begin{eqnarray}\begin{array}{} {{\tilde{\sigma }\,}^{(3);xxyy}} & = & {{\tilde{\sigma }\,}^{(3);yyxx}}\;,\quad {{\tilde{\sigma }\,}^{(3);xyxy}}={{\tilde{\sigma }\,}^{(3);yxyx}}\;,\quad {{\tilde{\sigma }\,}^{(3);xyyx}}={{\tilde{\sigma }\,}^{(3);yxxy}}, \\ {{\tilde{\sigma }\,}^{(3);xxxx}} & = & {{\tilde{\sigma }\,}^{(3);yyyy}}={{\tilde{\sigma }\,}^{(3);xxyy}}+{{\tilde{\sigma }\,}^{(3);xyxy}}+{{\tilde{\sigma }\,}^{(3);xyyx}}. \\ \end{array}\end{eqnarray} \tag{ 18 }$$

In nonlinear optics it is standard to symmetrize the terms ${{\tilde{\sigma }\,}^{(3);ijkl}}\left( {{\omega }_{j}},{{\omega }_{k}},{{\omega }_{l}} \right)$ by permuting the indices $\left( jkl \right)$ and the corresponding frequency variables $\left( {{\omega }_{j}},{{\omega }_{k}},{{\omega }_{l}} \right)$; [34] it is easy to confirm that such symmetrized coefficients, which we denote by ${{\sigma }^{(3);ijkl}}\left( {{\omega }_{j}},{{\omega }_{k}},{{\omega }_{l}} \right)$ can without impunity be used in place of ${{\tilde{\sigma }\,}^{(3);ijkl}}\left( {{\omega }_{j}},{{\omega }_{k}},{{\omega }_{l}} \right)$ in the expression for ${{J}^{(3)d}}\left( t \right)$, and they satisfy the symmetry relations given above.

The expressions given above are exact for the two-band model. However, for studying optical transitions around the Dirac points $\boldsymbol{k}=\frac{2\pi }{{{a}_{0}}}\left( \hat{\boldsymbol{x}}\,/\sqrt{3}+\hat{\boldsymbol{y}}\,/3 \right)$ or ${{\boldsymbol{k}}^{\prime }}=\frac{2\pi }{{{a}_{0}}}\left( \hat{\boldsymbol{x}}\,/\sqrt{3}-\hat{\boldsymbol{y}}\,/3 \right)$, it is usual to approximate the electronic dispersion relation as linear, i.e. ${{\epsilon }_{s\boldsymbol{k}+\boldsymbol{k}}}\approx s\hbar {{v}_{F}}k$ with Fermi velocity ${{v}_{F}}=\sqrt{3}{{a}_{0}}{{\gamma }_{0}}/2\hbar$, the velocity matrix elements as ${{\boldsymbol{v}}_{ss\boldsymbol{k}+\boldsymbol{k}}}\approx s{{v}_{F}}\;\boldsymbol{k}/k$ and ${{\boldsymbol{v}}_{s\bar{s}\boldsymbol{k}+\boldsymbol{k}}}\approx is{{v}_{F}}\hat{\boldsymbol{z}}\,\times \boldsymbol{k}/k$, and the interband Berry connection as ${{\xi }_{s\bar{s}\boldsymbol{k}+\boldsymbol{k}}}\approx \hat{\boldsymbol{z}}\,\times \boldsymbol{k}/2{{k}^{2}}$. These approximations follow immediately from equations (5), (6), and (9), and when they are made an analytic result is found for the linear optical conductivity,

$$\begin{eqnarray}&&{{\sigma }^{(1);xx}}\left( \omega \right)=\frac{i{{\sigma }_{0}}}{\pi }\left[ \frac{4|\mu |}{\hbar \omega }-G\left( \frac{\hbar \omega }{2|\mu |} \right) \right]\;,\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray*}&&G(x)=\ln \left| \frac{1+x}{1-x} \right|+i\pi \theta (|x|-1)\end{eqnarray*}$$

is a dimensionless complex function of a real variable x. Here $\theta \left( y \right)$ is the Heaviside step function, equal to 0 for $y<0$ and 1 for $y>0$. Compared to a full band structure calculation based on our tight-binding model, this analytic result gives an error less than 15% for $\hbar \omega <3$ eV. We also find an analytic result for the third order optical conductivity,

Here $c.p.$ stands for the simultaneous cycle permutation of indexes in $\left\{ {{\omega _1},{\omega _2},{\omega _3}} \right\}$ and $\left\{ {{A}^{(1)}},{{A}^{(2)}},{{A}^{(3)}} \right\}$ with

$$\begin{eqnarray*}&&{{A}^{(1)}}=\left[ \begin{array}{} -3 \\ 1 \\ 1 \\ \end{array} \right],\quad {{A}^{(2)}}=\left[ \begin{array}{} 1 \\ -3 \\ 1 \\ \end{array} \right],\quad {{A}^{(3)}}=\left[ \begin{array}{} 1 \\ 1 \\ -3 \\ \end{array} \right].\end{eqnarray*}$$

The linear conductivity is the same as found by, e.g. Gu et al [6], and reduces to the universal conductivity ${{\sigma }_{0}}$ when taking the Fermi energy $\mu \to 0$. We plot ${{\sigma }^{(1);xx}}$ in figure 1(a). Due to the gapless linear dispersion, the third order conductivities contain many possible divergences as the photon energies involved—including all the ${{\omega }_{i}}$ and their combinations—go to zero or to the doping induced gap ($2|\mu |$).

Zoom In Zoom Out Reset image size

In limit of zero doping, $\mu \to 0$ and $G\left( x \right)\to i\pi$; then ${{\sigma }^{(3);xxxx}}$ simplifies to equation (2), and ${{\sigma }^{(3);xxyy}}={{\sigma }^{(3);xyxy}}={{\sigma }^{(3);xyyx}}={{\sigma }^{(3);xxxx}}/3$. These terms all exhibit power-law divergences as photon energies vanish. For doped graphene, the chemical potential is only involved in the function $G\left( x \right)$, and a logarithmic divergence appears for suitable chemical potentials. This could be utilized as a flexible method of controlling the third order nonlinearities, by tuning the doping through adjusting a gate voltage or applying chemical doping. To present a summary of the results we scale all photon energies by the chemical potential and rewrite the third order conductivity as

$$\begin{eqnarray}&&{{\sigma }^{(3);dabc}}\left( \mu ;{{\omega }_{1}},{{\omega }_{2}},{{\omega }_{3}} \right)={{\sigma }_{0}}{{\left| \frac{1\text{eV}}{\mu } \right|}^{4}}{{S}^{dabc}}\left( \frac{\hbar {{\omega }_{1}}}{|\mu |},\frac{\hbar {{\omega }_{2}}}{|\mu |},\frac{\hbar {{\omega }_{3}}}{|\mu |} \right).\end{eqnarray} \tag{ 21 }$$

We show the coefficient ${{S}^{dabc}}\left( w,w,w \right)$ describing third harmonic generation in figure 1 (b), the coefficient ${{S}^{dabc}}\left( -w,w,w \right)$ describing two-photon absorption and Kerr effects for $w<2$ in figure 1(c), and typical components ${{S}^{dabc}}\left( -{{w}_{s}},{{w}_{p}},{{w}_{p}} \right)$ for parametric frequency conversion in figure 1(d). For the validity of the perturbation theory, the electric field should be less than a value ${{E}_{m}}$ at which the contributions to the optical current from the third order term are comparable to those from the linear term. Away from the divergent regime in figure 1, ${{S}^{dabc}}$ is usually of the order of ${{10}^{-19}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$, which gives an estimated ${{\sigma }^{(3)}}/{{\sigma }^{(1)}}\sim {{10}^{-15}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$ for a chemical potential $\mu \sim 0.1$ eV, which then leads to ${{E}_{m}}\sim 300\;\text{kV}\;\text{c}{{\text{m}}^{-1}}$. This value is at the high end of the range of electric field amplitudes used in experiments.

3.1. Third harmonic generation

For third harmonic generation we set ${{\omega }_{1}}={{\omega }_{2}}={{\omega }_{3}}=\omega$ and show the transitions at the left hand side of figure 2(a). From equation (20) the third order conductivity for this effect is given by

$$\begin{eqnarray}&&{{\sigma }^{(3);xxxx}}=\frac{i{{\sigma }_{0}}{{\left( \hbar {{v}_{F}}e \right)}^{2}}}{48\pi {{\left( \hbar \omega \right)}^{4}}}T\left( \frac{\hbar \omega }{2|\mu |} \right)\end{eqnarray} \tag{ 22 }$$

with $T\left( x \right)=17G\left( x \right)-64G\left( 2x \right)+45G\left( 3x \right)$; the other components are ${{\sigma }^{(3);xxyy}}={{\sigma }^{(3);xyxy}}={{\sigma }^{(3);xyyx}}={{\sigma }^{(3);xxxx}}/3$. If we take $\mu \to 0$, using (1) we find the effective bulk susceptibility is $\chi _{\text{eff}}^{(3);xxxx}=6.3\times {{10}^{-20}}{{\left( \frac{\text{eV}}{\hbar \omega } \right)}^{4}}{{\text{m}}^{2}}{{\text{V}}^{-2}}$. For nonzero $\mu$, ${{S}^{xxyy}}\left( w,w,w \right)$ is plotted in figure 1(b). The divergence at w = 2 is very weak compared to those at $w=2/3$ and w = 1, due to cancellations in T(x). Setting $\hbar \omega =0.72$ eV, as in the experiment by Kumar et al [15], our results are $\chi _{\text{eff}}^{(3);xxxx}=$ $3.25\times {{10}^{-19}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$ for $|\mu |\to 0$, and $6\times {{10}^{-19}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$ for $|\mu |=0.3$ eV ($w\approx 2.4$); these are both two orders of magnitude smaller than the experimental value. While our calculations do not include scattering and the effects of electron–electron interactions, there are a number of other possible reasons for this discrepancy, which include: (i) for third harmonic generation, the imaginary part of ${{\sigma }^{(3)}}$ has the same importance as the real part, and the contribution from $\boldsymbol{k}$ points far away from Dirac points may be not ignorable. (ii) Even around the Dirac points, the widely used linear dispersion approximation, which we have adopted here, may not be adequate for an accurate calculation of optical nonlinearities. Likely more important are the following reasons: (iii) for a very low doping level, the strong light intensity used in THG measurements may cause saturated absorption, especially if excited by high-repetition-rate laser pulses, this would lead to a nonlinearity of a totally different type, beyond the perturbation approach we apply here. (iv) Thermal effects may lead to a higher effective nonlinearity than the intrinsic one calculated here; their study is beyond the scope of the present work.

Zoom In Zoom Out Reset image size

3.2. Kerr effects and two-photon absorption

We now look at the frequency components ${{\omega }_{1}}=-{{\omega }_{2}}=-{{\omega }_{3}}=-\omega$. A direct calculation of ${{\sigma }^{(3);dabc}}\left( -\omega ,\omega ,\omega \right)$ leads to a divergent result which can be understood by rewriting equation (21) as

$$\begin{eqnarray}\begin{array}{} {{S}^{dabc}}(-w,w+{{\delta }_{1}},w+{{\delta }_{2}}) & = & \frac{S_{A}^{dabc}\left( w \right)}{{{\delta }_{1}}{{\delta }_{2}}}+\frac{S_{B}^{dabc}\left( w;{{\delta }_{1}} \right)}{{{\delta }_{2}}}+\frac{S_{B}^{dabc}\left( w;{{\delta }_{2}} \right)}{{{\delta }_{1}}} \\ {} & {} & +S_{C}^{dabc}\left( w;{{\delta }_{1}},{{\delta }_{2}} \right). \\ \end{array}\end{eqnarray} \tag{ 23 }$$

The divergent terms are related to the optical transition shown at the left hand side in figure 2 (b) and determined by

$$\begin{eqnarray}&&S_{A}^{dabc}\left( w \right)=-\frac{i{{\left( \hbar {{v}_{F}}e \right)}^{2}}}{12\pi {{w}^{2}}}H\left( \frac{w}{2} \right)\left( {{A}^{(1)}}+{{A}^{(2)}}+{{A}^{(3)}} \right)\end{eqnarray} \tag{ 24 }$$

and

$$\begin{eqnarray}&&S_{B}^{dabc}\left( w;\delta \right)=\frac{i{{\left( \hbar {{v}_{F}}e \right)}^{2}}\left( 3w+\delta \right)}{12\pi {{w}^{2}}\left( w+\delta \right)\left( 2w+\delta \right)}H\left( \frac{w}{2} \right)\left( {{A}^{(1)}}+{{A}^{(2)}}+{{A}^{(3)}} \right)\end{eqnarray} \tag{ 25 }$$

with $H\left( x \right)=G\left( x \right)+G(-x)=2i\pi \theta \left( |x|-1 \right)$. The term $S_{C}^{dabc}\left( w;{{\delta }_{1}},{{\delta }_{2}} \right)$ has no singularity as ${{\delta }_{1}},{{\delta }_{2}}\to 0$. Its components satisfy $S_{C}^{xyyx}\left( w;0,0 \right)=S_{C}^{xyxy}\left( w;0,0 \right)$ and

$$\begin{eqnarray*}&&\left[ \begin{array}{} S_{C}^{xxyy}(w;0,0) \\ S_{C}^{xyxy}(w;0,0) \\ \end{array} \right]=\frac{i{{\left( \hbar {{v}_{F}}e \right)}^{2}}}{24\pi {{w}^{4}}}\left\{ G(w)\left[ \begin{array}{} 32 \\ -32 \\ \end{array} \right]+K(\frac{w}{2}) \right\}\end{eqnarray*}$$

with

$$\begin{eqnarray*}&&K(x)=G(-x)\left[ \begin{array}{} 15 \\ -17 \\ \end{array} \right]+G(x)\left[ \begin{array}{} -61 \\ 35 \\ \end{array} \right]+\left[ 8x\frac{dG(x)}{dx}+2{{x}^{2}}\frac{{{d}^{2}}G(x)}{d{{x}^{2}}} \right]\left[ \begin{array}{} 1 \\ 1 \\ \end{array} \right].\end{eqnarray*}$$

We point out that the divergences as ${{\delta }_{i}}\to 0$ associated with ${{S}_{A}}$ and ${{S}_{B}}$ only appear for frequencies satisfying $\hbar \omega >2|\mu |$, where one-photon absorption is nonzero. So we focus on the frequencies satisfying $\hbar \omega <2|\mu |$, which is the regime where only two-photon absorption occurs, analogous to the situation in a gapped semiconductor with a gap corresponding to $2|\mu |$. In this regime, we find $\operatorname{Re}\left[ S_{C}^{xyxy} \right]=-\operatorname{Re}\left[ S_{C}^{xxyy} \right]$, and plot the w-dependence of ${{S}^{dabc}}(-w,w,w)=S_{C}^{dabc}\left( w;0,0 \right)$ in figure 1(c).

The two-photon absorption rate is written as

$$\begin{eqnarray}&&\frac{dn\left( t \right)}{dt}=\xi _{2}^{abcd}\left( \omega \right)E_{-\omega }^{a}E_{-\omega }^{b}E_{\omega }^{c}E_{\omega }^{d},\end{eqnarray} \tag{ 26 }$$

where n(t) is the areal density (for monolayer structures such as graphene) of electron-hole pairs injected. The coefficient $\xi _{2}^{abcd}(\omega )$ can be linked with the real part of ${{\sigma }^{(3);abcd}}(-\omega ,\omega ,\omega )$ by matching the energy loss of the light field $\left\langle \boldsymbol{J}(t)\cdot \boldsymbol{E}(t)\right\rangle$, written in terms of ${{\sigma }^{(3);abcd}}(-\omega ,\omega ,\omega )$, with the total energies of excited electron-hole pairs $\left\langle 2\hbar \omega \frac{dn\left( t \right)}{dt}\right\rangle$, written in terms of $\xi _{2}^{abcd}\left( \omega \right)$; we find $\xi _{2}^{abcd}\left( \omega \right)=3{{\left( \hbar \omega \right)}^{-1}}\operatorname{Re}\left[ {{\sigma }^{(3);abcd}}\left( -\omega ,\omega ,\omega \right) \right]$ yielding

$$\begin{eqnarray*}&&\xi _{2}^{xyxy}(\omega )=-\xi _{2}^{xxyy}(\omega )=\frac{4{{\sigma }_{0}}{{\left( \hbar {{v}_{F}}e \right)}^{2}}}{{{\left( \hbar \omega \right)}^{5}}}\theta (\hbar \omega -|\mu |),\text{for}\ \hbar \omega <2|\mu |.\end{eqnarray*}$$

This result agrees with calculations based on Fermi's Golden Rule [18] where only transitions at the right hand side of figure 2 (b) are considered and give $\xi _{2}^{dabc}\propto {{\left( \hbar \omega \right)}^{-5}}$ for $|\mu |<\hbar \omega <2|\mu |$.

In bulk materials the Kerr effect can be described by a nonlinear index of refraction ${{n}_{2}}$, leading to an index of refraction given by $n={{n}_{0}}+{{n}_{2}}I$ where I is the light intensity and ${{n}_{0}}$ is the index at very low intensities. With the very naïve model of graphene behaving as a thin layer of macroscopic material with effective response coefficients, we would have ${{n}_{2}}=3\chi _{\text{eff}}^{(3)}/\left[ 4{{\epsilon }_{0}}c\left( 1+\chi _{\text{eff}}^{(1)} \right) \right]$ [34]; in these formulas we assume that the light is linearly polarized, $\chi _{\text{eff}}^{(3);xxxx}$ is related to ${{\sigma }^{(3);xxxx}}(-\omega ,\omega ,\omega )$ by equation (1), and $\chi _{\text{eff}}^{(1);xx}$ to ${{\sigma }^{(1);xx}}\left( \omega \right)$ by the corresponding relation,

$$\begin{eqnarray*}&&\chi _{\text{eff}}^{(1);da}={{\sigma }^{(1);da}}/(-i\omega {{\epsilon }_{0}}{{d}_{gr}}).\end{eqnarray*}$$

Hence starting from our calculated quantities we can extract an effective ${{n}_{2}}$. In the experiment by Zhang et al [14], the photon energy was $\hbar \omega =0.8$ eV. Assuming $\mu \in \left[ 0.4,0.7 \right]$ eV, we calculate ${{n}_{2}}\sim {{10}^{-11}}\;\text{c}{{\text{m}}^{2}}\;{{\text{W}}^{-1}}$, much less than their measured value. However, if the doping level were much lower the experiment could have fallen in the divergent regime identified above, where we would not expect the perturbation approach adopted here would be sufficient.

3.3. Parametric frequency conversion

With a strong pump at frequency ${{\omega }_{p}}$, a signal photon at frequency ${{\omega }_{s}}$ can be converted to an idler photon at ${{\omega }_{i}}=2{{\omega }_{p}}-{{\omega }_{s}}$; this process is described by ${{\sigma }^{(3);dabc}}(-{{\omega }_{s}},{{\omega }_{p}},{{\omega }_{p}})$ and the corresponding optical transitions are illustrated at the right hand side of figure 2 (a). At $|\mu |\ll \hbar {{\omega }_{p,s,i}}$ we find the conductivity given by equation (2) is

$$\begin{eqnarray}&&{{\sigma }^{(3);dabc}}(-{{\omega }_{s}},{{\omega }_{p}},{{\omega }_{p}})={{\sigma }_{0}}{{\left( \hbar {{v}_{F}}e \right)}^{2}}/\left[ 6{{\hbar }^{4}}{{\omega }_{p}}{{\left( {{\omega }_{p}}-{{\omega }_{s}} \right)}^{2}}\left( 2{{\omega }_{p}}-{{\omega }_{s}} \right) \right].\end{eqnarray} \tag{ 27 }$$

This differs from the formula of Hendry et al [11] for the same process, although both this formula and that of Hendry et al exhibit an overall frequency dependence of ${{\omega }^{-4}}$. In the experiments of Hendry et al [11] the photon energies were $\hbar {{\omega }_{p}}=1.31\;\text{eV}$ and $\hbar {{\omega }_{s}}=1.02\;\text{eV}$ (corresponding to wavelengths of 950 nm and 1210 nm, and direct calculations from both formulas give similar values of ${{\sigma }^{(3);xxxx}}\sim {{10}^{-23}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$ (or $\chi _{\text{eff}}^{(3);xxxx}\sim {{10}^{-18}}\;{{\text{m}}^{2}}\;{{\text{V}}^{-2}}$), smaller by two orders of magnitude than the experimental values they extracted by comparing the nonlinear response with that of gold.

For large $|\mu |$, in figure 1(d) we plot the dependence of ${{S}^{xxyy}}\left( -{{w}_{s}},{{w}_{p}},{{w}_{p}} \right)$ and ${{S}^{xyxy}}\left( -{{w}_{s}},{{w}_{p}},{{w}_{p}} \right)$ on ${{w}_{s}}$ for a given ${{w}_{p}}=1.5$. For both components, divergences appear in the imaginary parts as ${{w}_{s}}\to 0$ and in the real parts as ${{w}_{s}}\to 2{{w}_{p}}$, and scale as $1/{{w}_{s}}$ and $1/\left( 2{{w}_{p}}-{{w}_{s}} \right)$ respectively. These divergences correspond to two interesting physical processes: as ${{w}_{s}}\to 0$, the idler frequency is close to twice ${{\omega }_{p}}$, which indicates that graphene could manifest very high second order nonlinearities in the presence of an external DC field [12]. For ${{w}_{s}}\to 2{{w}_{p}}$, the idler frequency vanishes, and the divergence ${{\left( 2{{\omega }_{p}}-{{\omega }_{s}} \right)}^{-1}}$ shows the existence of two-color coherent current injection, to which we turn in the next section. For ${{w}_{s}}\geqslant 1$, two photon absorption (or one photon absorption (${{w}_{s}}\geqslant 2$)) will occur and reduce the signal light intensity. However, the imaginary parts of ${{S}^{xyxy}}(-{{w}_{s}},{{w}_{p}},{{w}_{p}})$ show logarithmic divergences at ${{w}_{s}}\to 1$ and ${{w}_{s}}\to 2$, which could enhance the idler intensity.

3.4. Two-color coherent current injection

To link up with earlier treatments we now consider ${{\omega }_{1}}={{\omega }_{2}}=-\omega$ and ${{\omega }_{3}}=2\omega +\delta \omega$, as $\delta \omega \to 0$. The divergence, as ${{\left( \delta \omega \right)}^{-1}}$, indicates the injection of current through the interference of one- and two-photon absorption processes [35],

$$\begin{eqnarray*}&&\frac{d{{J}^{a}}}{dt}={{\eta }^{abcd}}(\omega )E_{-\omega }^{b}E_{-\omega }^{c}E_{2\omega }^{d}+c.c.\end{eqnarray*}$$

The tensor ${{\eta }^{abcd}}\left( \omega \right)$ can be calculated directly from Fermi's Golden Rule [18], but it also can be extracted from the expression for ${{\sigma }^{abcd}}\left( -\omega ,-\omega ,2\omega +\delta \omega \right)$,

$$\begin{eqnarray}&&{{\eta }^{abcd}}\left( \omega \right)=\underset{\delta \omega \to 0}{\mathop{\lim }}\,\left( -3i\left( \delta \omega \right){{\sigma }^{abcd}}\left( -\omega ,-\omega ,2\omega +\delta \omega \right) \right).\end{eqnarray} \tag{ 28 }$$

We find ${{\eta }^{xxyy}}={{\eta }^{xyxy}}$ and

$$\begin{eqnarray*}&&\left[ \begin{array}{} {{\eta }^{xxyy}}\left( \omega \right) \\ {{\eta }^{xyyx}}\left( \omega \right) \\ \end{array} \right]=\frac{i{{\sigma }_{0}}{{(\hbar {{v}_{F}}e)}^{2}}}{6\hbar {{\left( \hbar \omega \right)}^{3}}}\left\{ \theta (\hbar \omega -|\mu |)\theta (2|\mu |-\hbar \omega )\left[ \begin{array}{} 4 \\ -4 \\ \end{array} \right]+\theta (\hbar \omega -2|\mu |)\left[ \begin{array}{} 1 \\ 1 \\ \end{array} \right] \right\}.\end{eqnarray*}$$

Again, we can identify different regimes here: (i) $\hbar \omega <2|\mu |<2\hbar \omega$. Here only the first term contributes to the injection rates. This corresponds to the interference of one-photon absorption at $2\omega$ with two photon absorption at $\omega$ (see the transitions at the right hand side of figure 2 (c)), and corresponds to the usual current injection observed in semiconductors where $2\hbar \omega$ is greater than the band gap but $\hbar \omega$ is not. (ii) $\hbar \omega >2|\mu |$. In this regime one-photon absorption occurs at the fundamental frequency, and the second order contribution to current injection at the same position in the Brillouin zone leads to a divergent result. However, a new contribution to the coherent current injection appears and is finite; it arises from a new process of interfering pathways involve $\omega$ and $\omega -2\omega$ transitions [36] (at the left hand side of figure 2 (c)). For experiments in semiconductors, where the fundamental energy $\hbar \omega$ is typically less than the band gap, this regime is usually not explored. But note that in graphene the onset of this new channel leads to obvious amplitude changes in η, which could be easily detected in experiments by measuring the sudden change of the injection currents amplitude with scanning the fundamental frequencies across $2\mu$.

We investigated the linear and third order nonlinear optical conductivity of doped graphene, at zero temperature in the tight binding model, using the semiconductor Bloch equations. Analytic results were obtained around the Dirac points by using the widely accepted linear dispersion approximation. The third order optical conductivities exhibit a complicated dependence on the photon frequencies and chemical potential. We discussed in detail third harmonic generation, Kerr effects and two-photon absorption, parametric frequency conversion, and two-color coherent current injection. A nonvanishing chemical potential mimics many of the features that result from the presence of a band gap in normal semiconductors, and divergences can result when the energy of the effective gap is matched to any of the photon energies involved. The different third order processes considered exhibit a wide range of behavior, with each process having its own signature features. The important role played by the chemical potential allows for the generation of desired nonlinearities by electrically tuning or chemical doping; for low doping levels, the easily saturated absorption may lead to a nonlinear response totally different than those calculated here in the perturbative regime. Both systematic experiments and full band structure calculations, including finite temperature effects, thermal effects, scattering processes, and saturation, are clearly required to identify the true nature of the nonlinear response. But the perturbative calculations presented here identify a host of frequency regimes where very interesting behavior can be expected, even if more sophisticated calculations are required to elucidate it.

This work has been supported by the EU-FET grant GRAPHENICS (618086), by FWO-Vlaanderen which provides funding through the FWO project G.A002.13N and the Postdoctoraal Onderzoeker grant for N Vermeulen, by the Natural Sciences and Engineering Research Council of Canada, by VUB-Methusalem, VUB-OZR, and IAP-BELSPO under grant IAP P7-35.