{'fraction_non_alphanumeric': 0.04000487439722856, 'fraction_numerical': 0.01335236669394008, 'mean_word_length': 4.515506481036966, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Third order nonlinearity of graphene: effects of phenomenological relaxation and finite temperature

We investigate the effect of phenomenological relaxation parameters on the third order optical nonlinearity of doped graphene by perturbatively solving the semiconductor Bloch equation around the Dirac points. An analytic expression for the nonlinear conductivity at zero temperature is obtained under the linear dispersion approximation. With this analytic formula as starting point, we construct the conductivity at finite temperature and study the optical response to a laser pulse of finite duration. We illustrate the dependence of several nonlinear optical effects, such as third harmonic generation, Kerr effects and two photon absorption, parametric frequency conversion, and two color coherent current injection, on the relaxation parameters, temperature, and pulse duration. In the special case where one of the electric fields is taken as a dc field, we investigate the dc-current and dc-field induced second order nonlinearities, including dc-current induced second harmonic generation and difference frequency generation.


I. INTRODUCTION
The optical nonlinearities of monolayer graphene have recently attracted wide attention, 1-3 both experimentally and theoretically. The nonlinear susceptibility of graphene 4 is both strong -per atom it is orders of magnitude higher than that of common gapped semiconductors and metals -and controllable by the chemical potential, 5,6 which can be tuned by an external gate voltage 7,8 or chemical doping. 9 With the possibilities it offers for integration in silicon-based optical integrated circuits, graphene is an exciting new candidate for enhancing nonlinear optical functionalities in siliconbased on-chip optical devices, such as on-chip broadband light sources, electro-optic modulators, 10,11 optical switches, [12][13][14] and optical transistors. 15,16 In realizing some of these devices, 14 the presence of second order optical nonlinearities, especially second harmonic generation (SHG), is a key requirement.
Experimental studies of many of the optical nonlinear effects mentioned above have already demonstrated in graphene. Typically the experimental data are analyzed by extracting an effective optical nonlinear susceptibility, with the graphene monolayer treated as a thin film with a thickness of 3.3Å. [30][31][32] In this way, most of the experimental techniques used to determine the nonlinear optical response of bulk materials or thin films can be directly applied to the study of graphene. In a gapped semiconductor, third order susceptibilities do not change drastically in the nonresonant regime, where all photon energies are much lower than the energy gap. 33 Yet they show a strong and complicated photon energy dependence in pristine graphene because resonant transitions always exist for any photon energy, due to the vanishing gap and the presence of free carriers, leading to some similarities with a metal film. 32 These complexities have been observed in experimental studies of parametric frequency conversion, 30 THG, 31,32,34 Kerr effects and two photon absorption, 3,35-37 two color coherent control, [38][39][40] and SHG [17][18][19][20][21]27 in graphene.
Theoretically, the optical nonlinearities of graphene have been investigated by perturbative treatments based on Fermi's Golden Rule, and by density matrix calculations, both of which are standard methods in studying the optical response of gapped semiconductors. In an earlier communication we sketched some of the relevant work done before early 2014 41 ; recent contributions include a calculation by Mikhailov 42 of THG 43 . All of these studies focused on one or a few nonlinear effects. In our earlier work 5 we performed a perturbative calculation based on a density matrix formalism; ignoring all scattering effects, we obtained an analytic expression for the general optical sheet conductivity σ (3) (ω 1 , ω 2 , ω 3 ), which can be related to the effective susceptibility χ (3) (ω 1 , ω 2 , ω 3 ), in doped graphene at zero temperature. We found that the optical conductivities depend strongly on the chemical potential and photon energies, and exhibit many divergences associated with resonant transitions, which occur when photon energies or their combinations match the chemical potential gap. Taking ω 3 = 0 and including phenomenological relaxation times for the generation of both dc and optically induced current, we calculated the current induced second order nonlinearities at zero temperature and obtained an analytic expression 6 for CSHG. The effective susceptibility shows two peaks corresponding to two resonant transitions induced by the fundamental and the second harmonic light, with the peak values strongly dependent on the relaxation time. Adopting the parameters used in calculations of bilayer graphene, 28 we obtained a prediction of a peak susceptibility in monolayer graphene that was similar to that predicted for the bilayer; the EFISH contribution was ignored in that calculation.
The importance of the relaxation time demonstrated in that study, and the desire for more realistic calculations to compare with experiment, motivates the present work. Here we consider the inclusion of scattering effects in the semiconductor Bloch equations (SBE) within a relaxation time approximation, allow for finite temperature to the extent that it affects the initial state, and explicitly consider the nonlinear response to pulses of light. We obtain an analytic expression for the full nonlinear optical conductivity σ (3) (ω 1 , ω 2 , ω 3 ) for optical transitions around the Dirac points. We discuss the predictions that follow from this expression for different optical effects, and we compare with experiment where possible.
Our focus in this work is on doped graphene, where the chemical potential µ = 0. However, the chemical potential dependence of our general expression for σ (3) (ω 1 , ω 2 , ω 3 ) allows us to study the special case of µ → 0. At the very least we might expect that, for electrons close to the Dirac points, the distinction between "interband" and "intraband" motion could be lost. Although different terms that are nominally associated with interband and intraband motion arise naturally in the development of the perturbation series, the distinction between those two "kinds" of motion is at best approximate, 45 and we indeed find that the way those different formal terms contribute to the final result for small µ is nontrivial. More importantly, we generally associate the validity of a perturbative expansion of the optical response with the assumption that the energy induced by the presence of the optical field is much less than the energy difference between the bands. In graphene this is always violated for some states around the Dirac points, regardless of the strength of the optical field. If these states are occupied by electrons, as they are in undoped graphene, the reasonableness of a perturbative expansion is in doubt. Indeed, even a semiclassical treatment of the response to an applied electric field of electrons near the Dirac points exhibits a breakdown of the perturbative analysis 23 as µ → 0. We find evidence of the same kind of behavior in the quantum treatment presented here. This has consequences even for doped graphene if finite temperature is considered, for thermal fluctuations always place some electrons near the Dirac points.
We organize our paper as follows: In Section II we introduce the SBE and our approximations for including scattering effects; the details of the derivation of the nonlinear optical conductivity is given in Appendix A. The last two subsections of Section II address the extension of the calculation to finite temperature, and the treatment of the response to a pulse with finite duration. In Section III we discuss the third order nonlinear effects, including THG, Kerr effects and two photon absorption, two-color coherent current injection, and parametric frequency conversion; in Section IV we discuss the current-induced second order nonlinearities, including CSHG, EFISH, and the nonlinear optical conductivity σ (3) (−ω s , ω p , 0) that describes current-induced difference frequency generation. Throughout the sections we compare with experimental results when appropriate. We conclude in Section V.

II. MODEL
We take the Hamiltonian of graphene to be Here H 0 is the unperturbed electron Hamiltonian, where the a sk are annihilation operators of Bloch states |sk for band s and wave vector k, with eigen energy ε sk . Here H eR describes the interaction with radiation and in the dipole limit, where the electric field E(t) is approximated as uniform, we have where e = −|e| and is the Berry connection between states |s 1 k and |s 2 k , with A cell the unit cell area and u sk (r) the periodic part of the Bloch function, r|sk = (2π) −1 e ik·r u sk (r, z), where k = (k x , k y ) and r = (x, y); the graphene is assumed to lie in the x − y plane. We neglect any response of the system to electric field components in the z direction. The scattering terms are given by H ei for the electron-impurity scattering, H ep for the electron-phonon interaction, and H ee for the carrier-carrier scattering. The system is described by a density matrix that is initially diagonal both in band index and (continuous) where 0 ≤ n s1k1 ≤ 1 describes the initial occupation of the state. In the presence of an applied uniform electric field it remains diagonal in k 1 and k 2 but can acquire off-diagonal elements in s 1 and s 2 , describing the correlation between state amplitudes for |s 1 k 1 and . We can think of ρ s1s2k (t) as the elements of a 2×2 matrix ρ k (t), and their dynamics are determined by the SBE Here ∂ρ s 1 s 2 k ∂t scat includes the scattering terms induced by H ei +H ep +H ee , which could in principle be obtained from well-established treatments of many-particle systems, such as the many-particle density-matrix framework 46,47 or the Keldysh Green function method. 48 In an ordinary semiconductor with parabolic band structure, the current relaxation is mostly caused by carrier-phonon and carrier-impurity scattering, while carrier-carrier interactions are less significant due to the approximate equivalence of momentum conservation and velocity conservation. However, the novel linear band structure of graphene breaks this equivalence, and the carriercarrier interactions play an important role in current relaxation; 47,49 thus the full expression for the scattering terms is complicated and even hard to solve numerically. 48 We proceed in the standard way by assuming the validity of a perturbation expansion with ρ (n) s1s2k (t) = ρ 0 s1s2k = δ s1s2 n s1k is the density operator characterizing the equilibrium occupation of single-particle states at finite temperature T and chemical potential µ, where ρ (n) k ≡ 0 for n < 0. As a very rough approximation, a relaxation time approximation 49 can be adopted to give Here s1s2k is a relaxation parameter introduced to describe the dynamics of ρ s1s2k corresponds to a phenomenological relaxation time. In a real system, Γ (n) s1s2k can be expected to depend on the temperature, chemical potential, and external field. 50 Yet because the relaxation plays an important role in optical nonlinearities around resonant transitions, the extremely phenomenological treatment 51 in Eq. (8) can still reveal part of the physics, and in a very simple way. Even with the use of the six phenomenological constants Γ for interband transitions and Γ for intraband transitions, we are still able to obtain an analytic result for the perturbation calculation within the linear dispersion approximations around the Dirac points at zero temperature. From ρ (n) k (t), the (areal) current density, which in our model has only x and y components, is calculated as We give the derivation in Appendix A, where the spin degeneracy is included. We extract the linear optical conductivity σ (1);da (ω) from where E a (ω) = dtE a (t)e iωt . In graphene, the hexagonal lattice has D 6h (6/mmm) symmetry, 52 and there is only one independent nonzero component σ (1);xx = σ (1);yy . We first consider the zero temperature results.
In this paper, we restrict ourselves to the neighborhood of the Dirac points (see Fig. 1), assuming a linear dispersion relation with two relevant bands that we label s = + (upper) and − (lower). We recover the usual result 50,53 Here σ 0 = e 2 /(4 ) is the universal conductivity, and The third order current is given as Here the symmetrized third order optical conductivity σ (3);dabc is where the unsymmetrized third order optical conductivity is given as with e , ω 0 ≡ ω 2 + ω 3 , and ω ≡ ω 1 + ω 0 . We have followed the standard convention of nonlinear optics 33 in symmetrizing the terms σ (3);ijkl (ω j , ω k , ω l ) by permuting the indices (jkl) to arrive at the nonlinear conductivity σ (3);dabc (ω 1 , ω 2 , ω 3 ). The light-matter interaction in Eq. (3) can be formally separated into an interband contribution (s 1 = s 2 ) and an intraband contribution (s 1 = s 2 ), and the terms proportional to the different S i in σ (3);dabc can be classified according to how many times each contribution appears. 45 The term proportional to S 1 arises from only the intraband contributions, and the term proportional to S 6 arises from only the interband contributions; all others involve mixtures of both. The quantities S dabc i , σ (3);dabc , and σ (3);dabc are all fourth order tensors. Neglecting the optical response in the z direction, there are in all 8 nonzero components for the D 6h symmetry, among which three are independent; they are and In the following, we write the independent nonzero components of fourth rank tensors as column vectors, ordering the independent components of a fourth rank tensor where note A 0 = −(A 1 +A 2 +A 3 ), we can present the analytic expression for the different components of σ (3);dabc appearing in Eq. (15) at zero temperature, using the approximation of a linear dispersion relation around the Dirac points, as and where For the details see Appendix A.
All of these divergences only occur at complex frequencies in the presence of relaxation, and have their analogs in gapped systems. Of a different nature are the divergences that arise as |µ| → 0. While in a semiclassical calculation and in the absence of relaxation the intraband third order nonlinear response coefficient that can be extracted from the full nonlinear response is divergent 23 as |µ| −1 , one might hope that in the presence of relaxation this would be ameliorated. Yet in general it is not. To see this, we reorganize the unsymmetrized conductivity to write where σ A includes all terms involving G µ , σ B includes all terms involving H µ and I µ , and the remainder, σ C , includes all terms proportional to |µ| −1 . Similar separations are also used for the symmetrized conductivity σ (3);dabc . The term σ C can be simplified to yield Note that even for finite relaxation we havẽ σ (3);dabc C (ω 1 , ω 2 , ω 3 ) diverging as |µ| → 0, for general frequencies (ω 1 , ω 2 , ω 3 ), when Γ (j) e = Γ (j) i for both j = 1 and j = 3. At least within the simple description of relaxation we adopt here, the perturbation theory seems to demand that either the first or third order relaxation rates (or both) must not distinguish between intraband and interband relaxation to achieve a finite result as |µ| → 0. This is at least consistent with the physical intuition that the distinction between intraband and interband motion is blurred as |µ| → 0, in any case for electrons near the Fermi level, and any reasonable theory should respect that; recall that in our phenomenological description of relaxation all carriers share the same Γ i . But clearly a more sophisticated theory is in order to address the limit |µ| → 0.
More evidence for the blurring of the distinction between intraband and interband motion as |µ| → 0 can be seen from how the contributions to σ is the term proportional to S 1 ; it varies with |µ| as |µ| −1 , which is qualitatively different than the variation as |µ| of the corresponding Drude term in the linear conductivity. Yet as |µ| → 0 the contribution to σ (3);dabc (ω 1 , ω 2 , ω 3 ) involving only the formal interband (s 1 = s 2 ) component of Eq. (3), that is proportional to S 6 , also becomes important; while it includes contributions from G µ (θ 3 ) and G µ (θ), there is also a term proportional to |µ| −1 . The formally "mixed" terms, S j , with j = 2, 3, 4, 5, 7, 8, also provide terms proportional to |µ| −1 . The summation of all these terms, all formally involving different proportions of interband and intraband contributions, gives Eq. (31); the |µ| −1 behavior in σ C (ω 1 , ω 2 , ω 2 ) cannot be associated with motion that is just formally intraband. Now note thatσ C (ω 1 , ω 2 , ω 2 ) vanishes as all the Γ vanish. Yet here we physically would expect to recover the relaxation free, semiclassical result 23 of a perturbative response divergent as |µ| −1 , for ω i ≪ |µ|, and the result associated in that calculation with purely intraband motion. And we do recover it here, but in a nontrivial way: Although σ (3) C (ω 1 , ω 2 , ω 3 ) vanishes, when the other contributions to σ (3);dabc (ω 1 , ω 2 , ω 3 ) are assembled and the limit ω i ≪ |µ| taken we find in which the leading term is exactly the same as the contribution proportional to the S 1 term, and which agrees with the relaxation free, semiclassical calculation 23 involving only intraband motion. While the physically appropriate result of purely intraband, semiclassical motion is recovered in this limit as it should be, the connection to formally intraband, interband, and mixed responses in σ (3);dabc (ω 1 , ω 2 , ω 3 ) is less than direct. From Eqs. (30,31) we can also study more generally the limits as the relaxation rates are allowed to vanish. Here we discuss the simple case where the intraband and interband relaxation rates are the same for all orders, but perhaps different than each other: Γ (j) i = Γ i and Γ (j) e = Γ e . We find that as Γ i/e → 0 we recover from σ (3);dabc A (ω 1 , ω 2 , ω 3 ) the results derived earlier 5 in the absence of relaxation. We find that in this limit the contributions to σ (3);dabc B (ω 1 , ω 2 , ω 3 ) involving nonresonant transitions scale as Γ i . For resonant transitions, there are two cases that require further attention: (i) Taking ω comb to be a possible frequency combination appearing in the expression in Eqs. (19)(20)(21)(22)(23)(24)(25)(26), resonant transitions (real or virtual) occur as | ω comb | = 2|µ|. Then the function H µ or respectively, and then σ e ; its limit depends on the sequence of limits of Γ i → 0 and Γ e → 0, and so there seems to be no single well-defined relaxation free limit within this phenomenological theory. (ii) For some ω comb = 0, there can be divergences that occur at real frequencies in the absence of relaxation; below we discuss the behavior of σ (3);dabc near these divergences by considering the frequencies in the neighborhood of some of them.

B. Finite temperature
In calculating the response of a system to optical radiation, two effects of the temperature are usually considered: its role in establishing the initial electron distribution, and how it affects relaxation rates. In this work, the latter is implicit in our choice of relaxation rates. In our perturbative calculation, the former can be taken into account in the following simple way: Explicitly displaying the chemical potential and temperature dependence, we write n sk (µ, T ) for the electron distribution at equilibrium, and σ (3) (µ, T ) for the nonlinear conductivity. By using with F µ (x, T ) = [1 + e β(x−µ) ] −1 , the conductivity at finite temperature can be related to the zero temperature conductivity via Here the second line is obtained by using the partial integration and the condition σ ] is a pulse function located at x = µ with a width of the order of the thermal energy, the conductivity at finite temperature T can be obtained by averaging the zero temperature values over the chemical potential in an energy window with a width of the order of magnitude of the thermal energy. In a case where the chemical potential µ and the frequencies {ω i } are chosen to be away from resonant transitions, the conductivity is a smooth function around µ. Considering that the thermal energy k B T is only about ∼ 25.8 meV at room temperature, the conductivity at room temperature is close to the value at zero temperature away from resonant transitions. However, around resonant transitions where the conductivity diverges, the effects of finite temperature can be important. In Eqs. (19) to (26), the chemical potential appears in the functions G µ , H µ , and I µ in σ A and σ B , and as |µ| −1 in σ C . Therefore, the conductivity at finite temperature is determined by applying Eq. (34) to these quantities. The temperature effects on the contributions due to the functions G µ , H µ , and I µ are discussed in Appendix B.
Note that the treatment of the σ C term requires particular care, because at finite temperature there are always electrons initially near the Dirac points, and they will lead to the same prediction for divergent response that Eq. (31) indicates for electrons near the Dirac points at zero temperature in an undoped sample. To show this explicitly, from Eq. (34), |µ| −1 should be replaced by However, this diverges due to the singularity of the integrand at x = 0. Based on Eq. (31) where this term is nonzero only at Γ e , the divergence shows that either the perturbation theory or the assumption of unequal intraband and interband relaxation times in undoped graphene is not adequate, and more realistic treatments of the scattering and temperature are required. Nonetheless, from a full numerical solution of Eq. (5) and (8), 54 we find that contributions from the |µ| −1 term only give a small contribution to the total conductivity at finite temperature. Thus, at least at the level of the full SBE, whatever the final description of relaxation yields for the |µ| −1 term it will not lead to significant contributions. So for our finite temperature calculations we somewhat arbitrarily take C. Pulse response Because most nonlinear experiments are carried out using laser pulses, the optical response close to the divergences mentioned above is determined by the pulse shape. Except for the |µ| −1 divergences just discussed, the inclusion of the relaxation parameters Γ (n) i/e moves the divergent frequencies off the real axis. Yet it is necessary to investigate the pulse effects when the energy broadening of the pulse is larger than the broadening characterized by those relaxation parameters. For a field associated with pulses of a fixed polarization with the time domain envelope function p ωi (t), the Fourier transform is with the frequency domain envelope function The third order current in both time and frequency domain can be written as and We will be particularly interested in two special cases: (i) For δ i sufficiently small and σ (3);dabc sufficiently slowly varying in its frequency dependence so that over the frequency components of the envelope functions, the current response is given by In this case, the generated currents are also Gaussian in their time and frequency dependence.
(ii) For singular behavior where γ contains contributions from the relaxation parameters, the optical coefficient C dabc ω l ,ωm,ωn (t) satisfies The solution of this equation is and For a Gaussian pulse, we get , and Erf(x) is the error function. In the absence of relaxation, we have lim γ→0 This means that the current is nonzero even after the optical pulses have passed, indicating that current injection has occurred. For finite γ, C dabc ω l ,ωm,ωn (t) at t → ∞ is zero, but the injected current can still persist for some time. Fig. 2 shows the dependence of the current response on the pulse width. For a very long pulse, γ∆ ≫ 1, the current response has a shape that is nearly Gaussian; however, for γ∆ < 1, when the energy broadening of the pulse is larger than the relaxation rate, the current response obviously deviates from Gaussian shape, and can last long after the excitation pulses are passed.

III. THIRD ORDER OPTICAL NONLINEARITIES
To illustrate how relaxation affects the third order optical nonlinearities, in the sample calculations presented below we assume equal relaxation rates for all orders of response, putting Γ

A. Third harmonic generation
For monochromatic incident light with frequency ω, light is nonlinearly generated to lowest order at the third harmonic frequency 3ω and at the fundamental frequency ω. The first is described by the conductivity σ (3);dabc (ω, ω, ω); the second corresponds to Kerr effects and two photon absorption, both described by σ (3);dabc (−ω, ω, ω), and can be considered as a nonlinear correction to the linear optical response. In this section, we consider THG.
The conductivity tensor for THG only has one independent component The induced current responsible for the THG is Fig. (3) we give the result for σ (3);xxxx (ω, ω, ω) at |µ| = 0.3 eV for zero and room temperature. We first look at the results for zero temperature. The relaxationfree results are given as the thick (red) curves in Fig. 3 (a). This figure shows the step function of the real parts and the logarithmic divergence of the imaginary parts at three resonant photon energies ω = 0.2, 0.3, and 0.6 eV, which correspond to the resonant transitions for which the chemical potential gap 2|µ| matches the energies of three photons, two photons, and one photon, respectively. 5 With relaxation included, the conductivity is a smooth function of ω, and plotted in Fig. 3 (b), (c), and (d). Some common effects induced by the relaxations are shown: (i) the divergent peaks of the imaginary parts of the conductivity in Fig. 3 (a) become finite and broadened, (ii) the step functions of the real parts become continuous, (iii) the real parts become finite as ω < 2|µ|/3 = 0.2 eV, and increase rapidly with decreasing frequency. They receive contributions not only   from intraband transitions, describing Drude-like effects, but also from the interband transitions due to the linear dispersion relation of graphene (for example, see the prefactor ( ω) −4 in Eq. (52)).
To illustrate the dominant features in these fine structures, we can analytically expand the coefficients of the functions G µ , H µ , I µ , and the |µ| −1 term in the conductivity, for small relaxation parameters Γ i,e /( ω) ≪ 1, to write In the relaxation-free limit as Γ e,i → 0, σ ;xxyy B → 0 and σ (3);xxyy A recovers the results of our previous work. 5 However, the relaxation-free limit of σ (3);xxyy C strongly depends on the details of the chemical potential and relaxation parameters; this is the contribution to σ (3);xxyy C from the general term discussed earlier in Eq. (31), which is problematic unless Γ i = Γ e . For doped graphene where µ is finite, σ (3);xxyy C goes to zero with decreasing relaxation parameters; for graphene that is undoped or at low doping, a more sophisticated treatment is in order, as discussed in Section II A. For the limit Γ i,e , |µ| ≪ ω, the THG coefficient is approximated as The term proportional to |µ| −1 did not arise in our previous calculation 5 , where we assumed that Γ i,e → 0 faster than µ → 0. Deferring the treatment of small doping to later studies, we focus here on graphene with large enough chemical potential that σ (3);xxyy C (ω) does not make a significant contribution to the full third harmonic conductivity.
At room temperature, the conductivities for different relaxation parameters look very similar to each other, and the fine structures caused by the resonant transitions are smeared out. This can be understood by the results in Appendix B: temperature affects the conductivity by smearing and lowering the peaks caused by functions G, H, and I, which has an effect similar to increasing the value of Γ e . If we increase each Γ e by the thermal energy of room temperature, the values of these new Γ e in the four cases presented in Fig. 3 are close, and it is not surprising that we get similar room temperature results.

B. Kerr effects and two photon absorption
We now turn to the light nonlinearity generated at the same frequency ω of the incident light.
and consider The nonlinear response at frequency ω is then given by For linearly polarized light (θ = 0), the second term vanishes; the current from the first term has the same polarization as the incident field, and gives an intensity dependent correction of the linear conductivity σ xx eff (ω) = σ (1);xx (ω) + σ nl (ω), with An effective nonlinear susceptibility can be introduced 30-32 χ nl (ω) = σ nl (ω)/(−iωǫ 0 d gr ), where the effective thickness of graphene single layer d gr is taken to be 3.3Å; 30 from this an effective nonlinear refractive index n 2 and nonlinear loss β TPA can be extracted. In general, σ (1);xx has both real part and imaginary parts, and the calculation of n 2 and β TPA should follow the results of del Corso and Soles [55].
In the limit of no relaxation, we showed earlier 5 that σ (3);dabc (−ω, ω, ω) has many divergences, and its behavior in the neighborhood of equal frequencies can be written as where T dabc , which are real functions and are only nonzero when the photon energy is greater than that for the onset of one-photon absorption ( ω > 2|µ|). For fixed photon energy ω, the appearance of the divergence as µ decreases from 2|µ| > ω to 2|µ| ≤ ω indicates that it is associated with the existence of electrons (holes) at the k where one-photon absorption is possible. Physically, at these k the third order correction to one-photon absorption would lead to the perturbative description of the saturation, but in the absence of relaxation that correction diverges, as it would for an inhomogeneously broadened collection of two-level systems. At zero temperature, the sharp Fermi surface can strictly exclude such electrons (holes) for 2|µ| > ω. However, at finite temperature thermal fluctuations will always place some electrons (holes) at k where one-photon absorption can occur, and so the divergence in the third order response will exist for any photon energy.
(2) |µ| < ω < 2|µ|, in which two-photon absorption is present but one-photon absorption is still absent. The real part of the conductivity here scales as ( ω) −4 . Around ω = 2|µ|, the imaginary part shows a divergence ( ω − 2|µ|) −2 . For frequencies satisfying ω < 2|µ|, if the graphene is subject to a Gaussian pulse sufficiently narrow in frequency, the nonlinear current induced will still have a shape that is approximately Gaussian, and characterizing the nonlinear response to a pulse by Eq. (56) makes sense. (3) ω > 2|µ|, where both twoand one-photon absorption are present. The imaginary part of the conductivity diverges as ( ω − 2|µ|) −2 around ω = 2|µ|, and the real part diverges for the entire region ω > 2|µ|. At finite temperature and in the absence of relaxation, the real part diverges for any photon energy ω. As we discussed after Eq. (57), the divergence of the real part of σ (3);dabc (−ω, ω, ω) is induced by the existence of electrons (holes) at the k where one-photon absorption occurs; at zero temperature, these electrons (holes) only exist when the chemical potential |µ| < ω/2, while at finite temperature, they exist at any chemical potential due to thermal fluctuations.
In Fig. 4 (b), (c), and (d), we present the results for the same relaxation parameters as those adopted in the THG calculation. Relaxation affects the conductivity in a complex way, but there are some qualitative features that can be identified: (i) In the neighborhood of the divergences that arise in the relaxation-free calculation, including the divergent regime ω > 2 |µ| and the special frequency ω = |µ|, both the real and imaginary parts of the conductivity are lowered and are everywhere finite. In Fig. 4 (b), (c), and (d), we find that a larger Γ e gives lower and broader peaks at ω = |µ| and 2|µ|. (ii) For the relaxation parameters used here, the real part of the conductivity is negative for ω > 2|µ|. Because of the presence in this frequency range of one-photon absorption, which is always positive, the two-photon absorption processes indicated by the real part of σ (3);dabc (−ω, ω, ω) can be understood as a correction to the simple linear prediction of the absorption. In fact, we can find a range of electric fields large enough so that σ da eff (ω) is negative; for a field anywhere near or above this strength the perturbative result is naturally suspect. (iii) Even for frequencies in the range |µ| < ω < 2 |µ|, where only two-photon absorption is present in the absence of relaxation, the real part of the nonlinear conductivity σ (3);dabc (−ω, ω, ω) can be negative. Yet in the presence of relaxation the linear conductivity σ (1);xx (ω) acquires a real part in this frequency range, and the real part of σ xx eff (ω), for example, is always positive for small enough electric field amplitudes, indicating absorption. However, these results indicate the sensitivity to the relaxation parameters of both the third order conductivity, and its interplay with the first order conductivity, and a more sophisticated description of the scattering is clearly in order.
At room temperature, the peaks or divergences are further broadened. For a given frequency ω, the regime ω > 2|µ| always contributes to a finite temperature calculation due to the average over the chemical potential.
The absolute values of the real part of the conductivity in the regime ω < 2|µ| also significantly increase.

C. Two-color coherent current injection
Now we turn to two-color coherent current injection, with frequencies ω 1 = ω 2 = −ω and ω 3 = 2ω. In the relaxation-free calculation, the conductivity σ (3);dabc (−ω, −ω, 2ω) diverges, and it is the divergence that describes the current injection. In fact, in the neighborhood of these frequencies, the conductivity can be written as where the injected current is determined by a wellbehaved function η dabc (ω), and σ (ω) is a smooth function of δω. With the inclusion of relaxation, the conductivity σ (3);dabc (−ω, −ω, 2ω) itself is well behaved. The divergence term in the relaxation-free limit becomes a term similar to the right hand side of Eq. (45). To check whether Γ i and Γ e have the same importance for the injection, we give the pulse calculations of C xxxx −ω,−ω,2ω (t) (see Eq. 39) in Fig. 5 for different relaxation parameters. After the laser pulse, the current response persists for  times associated with Γ i , showing that the contribution from Γ i dominates the relaxation of the injected current, as might be expected. To highlight this, we write From Eq. (15), only terms S dabc 1−4 including (δ 1 + δ 2 + δ 3 + iΓ i ) −1 contribute to η dabc (ω). By writing σ dabc inj (ω) ≡ η dabc (ω)/1eV, the first term in Eq. (64) becomes  where terms proportional to |µ| −1 are neglected. As Γ i , Γ e → 0, only terms involving Im[G µ (ω)] ∝ θ(|ω|−2|µ|) remain; σ dabc inj (ω) is a pure imaginary quantity, and is consistent with our previous work. 5 With relaxation included, terms involving Re[G µ ] and H µ appear.
In Fig. 6 (a) and (b) we plot the photon energy dependence of σ xxxx inj (ω) for different relaxation parameters at zero and room temperature. The real part is much smaller than the imaginary part. At zero temperature, we see in Fig. 6 (b) that for Γ e = 0.5 meV there are fine structures in the spectrum of Im[σ xxxx inj (ω)] around ω = |µ| and ω = 2|µ|; it is due to the H µ terms. As can be seen in the inset of Fig. 6 (b), finite temperature and finite Γ e lead to similar broadening and lowering of the peaks.
We also plot σ xxxx R (ω) in Fig. 6 (c) and (d). The amplitude of σ xxxx R is of the same order of magnitude as that of σ xxxx inj . However, due to the prefactor 1eV/(3Γ i ), which relates σ xxxx inj to its contribution to σ (3);dabc (−ω, −ω, 2ω) in Eq. (65), σ xxxx inj dominates for small Γ i , usually taken to be a few tens of meV.

D. Parametric frequency conversion
Third order nonlinearities can lead to the appearance of new frequencies via parametric frequency conversion, which is described by σ (3);dabc (−ω s , ω p , ω p ). Here ω p is the frequency of a strong pump field, and ω s is the signal frequency converted by interaction with the pump to an idler frequency ω i = 2ω p − ω s . Possible resonant transitions occur as any of the frequencies |ω p |, |ω s |, |ω p − ω s |, |2ω p |, or |2ω p − ω s | equal 2|µ|/ . In Fig. 7 we plot the dependence of σ (3);xxxx (−ω s , ω p , ω p ) on ω s for different relaxation parameters at ω p = 0.8 eV and |µ| = 0.5 eV.
(2) Around ω s = ω s2 + δω, the signal frequency is close to the onset of the one-photon absorption. For non-resonant transitions in a usual semiconductor, ω s2 and ω s1 are interchangeable frequencies to give the same conductivity of parametric frequency conversion; 33 here in graphene they yield asymmetric peaks because the resonant transitions dominate. In the limit of no relaxation, the conductivity shows a logarithmic divergence that is easily smeared out by the inclusion of small relaxation parameters.
(4) At finite temperature there is a further smearing of the peaks around the resonances, as described in Appendix B.
Besides these resonant transitions, two singularities are apparent: (i) the singularity around ω s = 2ω p , which corresponds to two-color coherent current injection. The singularity is not determined by the behavior of the functions G, H, I, but by the coefficients that premultiply them in Eqs. (19) to (26). Around ω s = 2ω p we put ω s = 2ω p + δω and find an equation similar to Eq. (64) discussed in Section III C. Two-color coherent current injection requires both onephoton absorption (for ω s ) and two photon absorption (for ω p ), i.e., ω p > |µ|. The parameters we have adopted in Fig. 7 fulfill this criterion, and thus the singularity appears. Finite temperatures do not qualitatively affect this singularity because it is not related to the chemical potential. (ii) The strong response around ω s = ω p is related to the third order correction to one-photon absorption, and only appears at finite temperature. Since |µ| < ω p < 2|µ|, at zero temperature only two-photon absorption is present and there is no one-photon absorption. However, at finite temperature thermal fluctuations will place electrons where one-photon absorption can occur, and the third order correction to that will lead, in the absence of relaxation, to a divergent result as discussed in Section III B; in the presence of relaxation the result will not be divergent but very large, describing the saturation of the one-photon absorption at the level of the third order response.

E. Comparison between calculations and experiments
Experiments have already extracted values of the effective third order susceptibilities of THG, 31,32,34 twophoton absorption, 34,36 Kerr effects, 3,35,37 and parameter frequency conversion 30 at some photon energies. The nonlinear conductivities we have calculated here are related to the effective susceptibility by 5,30 We first look at the THG, for which the experimental technique is perhaps the most mature, and the extracted values can likely be considered more reliable than those from other effects. For a reasonable chemical potential estimated from the sample preparation, the calculations without relaxation parameters 5 yield theoretical results for the nonlinear conductivity about two orders of magnitude smaller than the value extracted from experiments.
Here we have found that calculations at finite temperature for different sets of relaxation parameters (see the insets of Fig. 3) are almost the same as calculations at zero temperature and neglecting relaxation. For Kerr effects, because of the existence of divergent terms in the expressions and the probably very low chemical potential in experiments, it is not surprising that we could fit the nonlinear susceptibility at one photon energy by tuning the relaxation parameters. The complicated dependence is shown in Fig. 4. As ω > 2|µ|, the nonlinear conductivity at both zero and room temperatures can vary many orders of magnitude, depending on the relaxation parameters adopted.
For parametric frequency conversion observed in the experiment by Hendry et al., 30 with parameters ω p = 1.31 eV, ω s = 1.05 eV, and assuming a low chemical potential |µ| = 0.1 eV, we checked the dependence of the conductivity on the relaxation parameters Γ i and Γ e in the range of [0, 60] meV. We find the dependence is weak and the calculated values are still smaller than their claimed values by two orders of magnitude. 30 Admittedly the measured effective susceptibilities for parametric frequency conversion, Kerr effects and two photon absorption, and THG show a strong dependence on the measurement method, light frequency, pulse duration, and perhaps sample preparation. Yet even taking this into account, the conclusion that the theoretical results are about two orders of magnitude smaller than the measured results is inescapable. These discrepancies could arise for a number of reasons, including: (1) The samples in many experiments are not suspended graphene, but graphene on a substrate or in solution. Thus there may have been contributions to the optical nonlinearity from the interaction between the graphene sheet and its environment, which may be crucial considering that graphene is a one-atom thick material. (2) Thermal effects 33 caused by a high repetition rate of laser pulses, as used in Z-scan experiments, may play an important role. 56,57 (3) Because of the zero gap of graphene and the intense laser beams used in experiments, saturation 58 induced by one and/or two photon absorption can make necessary a treatment more sophisticated than that of perturbation theory. Zhang et al. 59 used the density matrix method to study four wave mixing in undoped graphene in the saturation regime, and found an effective χ (3) eff about 10 −17 m 2 /V 2 , and decreasing with increasing light intensity. Additional calculations for different third-order nonlinear effects in graphene in the saturation regime are needed to assess the impact of saturation on the theoretical nonlinearities. (4) The calculation at the independent particle level, which works well as a starting point for most gapped semiconductors, may fail in graphene, and it may be necessary to do a more realistic calculation, including the full band structure, and the detailed effects of scattering and the electron-electron interactions.

IV. DC CURRENT INDUCED SECOND ORDER NONLINEARITY
We now turn to the limiting case where one of the electric fields is a dc field, taking ω 3 = 0. The calculation of σ (3);dabc (ω 1 , ω 2 , 0) from Eq. (15) includes a term proportional to Therefore a nonzero relaxation Γ for the dc field is necessary to set up a steady state with a dc charge current in graphene. For other transitions included in σ (3);dacb (ω 1 , 0, ω 2 ) and σ (3);dcab (0, ω 1 , ω 2 ), it is not necessary to include relaxation associated with the dc field, because the dc field acts on the optical excitation with frequency ω 2 and ω 1 + ω 2 respectively; these only survive during the optical pulse. We list the relaxation parameters used in calculating the unsymmetrized conductivities in Table I. The third order conductivity of interest here, which we can refer to as the dc-induced second order optical conductivity, can be written as unsymmetrized σ relaxation parameters unsymmetrized σ relaxation parameters Γ  The first term includes all contributions that diverge as 1/Γ dc i , which are only involved in calculating σ (3);dabc (ω 1 , ω 2 , 0) and σ (3);dbac (ω 2 , ω 1 , 0); they both occur with the dc charge current. Thus we can associate it with the dc current-induced second order conductivity, in that it is second order in the optical fields at ω 1 and ω 2 . The second term includes all other contributions, and we can associate it with a dc field-induced second order conductivity, which exists even for a gapped semi-conductor without doping. Examining Eq. (15), we see that σ dabc J (ω 1 , ω 2 ) is independent of Γ dc i and Γ dc e , and it can be written as As we show below, the values of σ J and σ E are typically of the same order of magnitude; hence it is the value of Γ dc i that determines whether the dc-current induced second order conductivity or the dc-field induced second order conductivity makes the larger contribution to the dc-induced second order conductivity σ (3)dabc (ω 1 , ω 2 , 0). We can get a rough estimation of Γ dc i from the graphene mobility µ m . The dc limit of the optical conductivity can be obtained from Eq. (11) as The connection between the mobility and conductivity can be written as σ (1);xx (0) = N e |e|µ m with the carrier density N e = |µ| 2 π( vF ) 2 obtained from the linear dispersion.
Hence, we get For a sample with mobility µ m = 10 3 cm 2 /(V·s) and chemical potential µ = 0.5 eV, Γ dc i is about 10 meV. We will see below that for samples with such mobilities the dc-current induced effects will typically dominate the dc induced second order conductivity.
From Eq. (70), we see that the first resonance in σ dabc J (ω) is determined by H µ (2 ω + iΓ op e ) and I µ (2 ω + iΓ op e ); the other resonance is determined only by H µ ( ω+ iΓ op e ). Obviously, for both transitions smaller values of Γ op e result in a larger value and a sharper peak. At room temperature, these peaks are broadened and lowered. The vertical line at ω = |µ| in Fig. 8 (a1) comes from Im[H µ (2 ω + iΓ op e )] as Γ op e → 0, which is proportional to δ( ω − |µ|); the peak in Fig. 8 (c1) shows the fine structure of Im[I µ (2 ω + iΓ op e )] for Γ op e = 0.5 meV (see Appendix B). However, it is interesting to note that these two fine structures undergo important changes at room temperature: In Fig. 8 (a1), we see that the first fine structure leads to a peak with broadened width; in Fig. 8 (c1), we see that the second fine structure leads to a sign change around ω = |µ| when the temperature increases from zero to room temperature.

B. DC current induced difference frequency
A counterpart of the third order parametric frequency conversion discussed in Section III D is difference frequency generation which is, like second harmonic generation, a second order nonlinear effect that can be induced in graphene when applying a dc field. With a strong pump at frequency ω p , difference frequency generation converts a signal frequency ω s to a new frequency ω p −ω s ; the response is determined by σ (3);dabc (−ω s , ω p , 0). Similar to dc-induced second harmonic generation, there are current and electric field contributions to dc-induced difference frequency generation. As we found in Section IV A for dc-induced second harmonic generation, the current contribution should dominate the dc-induced difference frequency generation in a high mobility sample. As an example, we plot the chemical potential dependence of σ (3);dabc J (−ω s , ω p ) for different optical relaxation parameters in Fig. 9 for ω p = 1 eV (with a wavelength of about 1.24 µm) and ω s = 0.58 eV (with a wavelength of about 2.1 µm). For vanishing optical relaxation parameters (Γ op i = Γ op e = 0), it is clear from Fig. 9 (a) that there are three resonant transitions in the plotted chemical potential range: |µ 1 | = ( ω p − ω s )/2 = 0.21 eV, |µ 2 | = ω s /2 = 0.29 eV, and |µ 3 | = ω p /2 = 0.5 eV.
Without relaxation, the imaginary part of the conductivity is always zero except at these three resonant transitions (shown as vertical lines): the first is given by lim δ→0 , the other two are given by lim With finite relaxation rates or at finite temperature, the vertical lines are broadened to structures of finite strength and width.

C. Comparison between calculations and experiments
Bykov et al. 19 observed that SHG radiation from a graphene/SiO 2 /Si(001) substrate strongly depends on the applied current density in the graphene layer, which is attributed to the CSHG effect of graphene. A similar structure was also studied by An et al., 20,21 who could measure the radiation from different locations on the graphene sheet; they interpreted the result as EFISH, where the electric field is induced by current-associated trapped charge at the graphene/SiO 2 interface. Because of the interface contribution to the SHG radiation, [17][18][19][20][21] the contribution of the current related SHG from the graphene is hard to extract.
The best way of measuring the dc-induced second order nonlinearity of graphene, without any background contribution from interface effects, would be to mount graphene in a symmetric structure; this can be difficult. However, within the framework of the experiments of the type that have already been done, we can suggest a strategy that might help identify the in-plane graphene CSHG(EFISH) by the azimuthal angle dependence of the generated signal. For linearly polarized light with We see that the Cartesian components of the induced current vary as cosinusoidal functions of 0φ (that is, independent of φ) and 2φ. In a short-hand notation, we will characterize these as 0φ and 2φ dependences. In most experiments, 17-21 the graphene sample is mounted on a SiO 2 /Si substrate, where the interface between SiO 2 and Si gives an interface-induced SHG and the bulk Si gives an electric quadrupole/magnetic dipole induced SHG. But for different crystal orientations of the Si substrate, the dependence of the combined interface and bulk contributions on azimuthal angle will be different: 22 For the (111) face, the second harmonic radiation depends on the angle as 0φ and 3φ; for the (001) face, the dependence is 0φ and 4φ; while for the (110) face, the dependence becomes 0φ, 2φ, and 4φ. Therefore, from the azimuthal dependence of the SHG signal it might be possible to distinguish the graphene CSHG (EFISH) from the interface contributions, for example by putting graphene on top of different SiO 2 /Si structures, one with the (111) face of Si normal to the interface and one with the (001) face normal. Because of the same origin of CSHG and EFISH in graphene, they would have the same angle dependence, so such experiments would not help to distinguish between these different contributions from graphene; but for a heavily doped and high mobility graphene sample our calculations show that the CSHG should dominate.

V. CONCLUSION AND DISCUSSION
Perturbative analyses play a central role in nonlinear optics. Even when the electrons in a material are treated as independent, and relaxation is only described phenomenologically, the calculated response tensors that relate the induced polarization or current to powers of the applied fields indicate the nonlinear optical effects that are allowed, and point to where resonances can lead to an interesting dependence on time and frequency. Often more sophisticated models of the electron dynamics are required, and sometimes the perturbative framework itself is insufficient to address the physics of interest. But even then these kinds of perturbative treatments provide a starting point for more realistic calculations.
In this paper we have provided such a treatment of the nonlinear third order optical response of doped graphene, with the main goal of investigating the effects of phenomenological relaxation parameters, finite temperature, and laser pulse width on the induced currents. We focused on the contributions of optical transitions around the Dirac points, where the widely used linear dispersion relation is a good approximation. By solving semiconductor Bloch equations perturbatively, an analytic expression for general third order conductivities was obtained at zero temperature, taking different relaxation parameters for interband and intraband optical transitions. The nonlinear conductivities at finite temperature were obtained by an appropriate integration over the chemical potential. The conductivities show a complicated dependence on photon energy, chemical potential, and the relaxation parameters.
Even with the inclusion of relaxation we found that the perturbative approach itself is problematic at vanishing chemical potential, as might be expected from a similar result in the semiclassical limit, 23 except in the special case that either first-or third-order interband and intraband relaxation rates are set equal. The perturbative approach adopted is unproblematic for doped graphene at zero temperature, but is a concern at finite temperature, since thermal fluctuations always place some electrons or holes near the Dirac points. Yet numerical calculations of the full semiconductor Bloch equations indicate that the contribution of such electrons to the full optical response is small, so this effect does not afflict our results.
We discussed in detail different nonlinear effects, including third harmonic generation, Kerr effects and twophoton absorption, two-color coherent current injection, parametric frequency conversion, and dc-current andfield induced second harmonic generation and difference frequency generation. The interband relaxation generally broadens and lowers the resonant peaks, while the intraband relaxation plays an important role in some of the effects, including two-color coherent current injection and the dc-current induced second order nonlinearities. At room temperature most of the resonant structures are smeared out.
We also considered the response of graphene to laser pulses. The optical response depends in detail on the frequency width of the incident pulse and the frequency structure of the response tensors. The two natural limits are (1) when the frequency structure of the response coefficients is rather flat on the scale of the frequency width of the incident pulse, and the induced current follows the injecting pulses, and (2) when there are divergences in the response coefficients at frequencies close to the real axis, as for two-color coherent control, the a dynamics associated directly with the relaxation processes.
Comparison of our results with experiments is difficult, since in many of the reported experiments the graphene samples have not been characterized in the linear regime, and neither the relaxation parameters nor even the chemical potential have been identified. Results for some nonlinear response coefficients, such as that describing the Kerr effect, are predicted by our calculations to be so sensitive to these parameters that we cannot hazard a comparison of theory to experiment. Yet our results for third harmonic generation and parametric frequency conversion are insensitive enough to these parameters that we can conclude our results are about two orders of magnitude smaller than those extracted from experiments, 5 even with the adoption of reasonable relaxation parameters. We speculated on the causes of this disagreement; it is of course early days for both detailed experimental and theoretical studies of such nonlinear effects in graphene. But these disagreements may persist, and a long and difficult journey may be necessary to understand the details of the full nonlinear optical response of graphene. Even so we can expect that, as in the study of the nonlinear optical response of other materials, the kind of perturbative calculation we have presented here will provide a useful port of embarkation, paving the way to a better physical insight in the complex nonlinear optical response of graphene.
We describe the electronic states in graphene by the tight binding model employing carbon 2p z orbitals with only nearest neighbor coupling. Neglecting the overlap between different p z orbitals, the band structure that results is electron-hole symmetric with energies ε +k = −ε −k where +(−) is the band index for π * (π) bands, and the Berry connections satisfy ξ ssk = ξs sk ≡ r k and ξ ssk = ξss k where s = ± is the band index ands is the index of the band that is not the s band. 5 Up to the dipole approximation of light-matter interaction, the SBE in Eq. (7) can be expanded as with ǫ k = ε +k − ε −k . In Eq. (A1) the terms involving r k give the interband contribution, and the terms involving ∇ k give the intraband contribution. Treating the electric field term perturbatively, the first three terms are expanded as ρ with ω 0 = ω 2 + ω 3 and ω = ω 1 + ω 0 . By substituting the above expansion into Eq. (A1), we get the following equations for P (i) .
(1) The linear order terms are determined by Here v s1s2k are the matrix elements of the velocity operator, which satisfy v ++k = −v −−k and v +−k = −v −+k in the tight binding model we have adopted.
(2) The second order terms are determined by ν 0 P The solutions are P (2);bc Because the graphene crystal structure is centrosymmetric, its second order conductivity is zero.
The frequency dependence of P Here terms P abc i;sk with i = 1, · · · , 4 are related to the populations at band s and are given by to expand the derivatives in S i , and find that all the required integrations over k can be related to Using partial fractions and taking the integral over Ω, we obtain all the S i given in Eq. (19)(20)(21)(22)(23)(24)(25)(26)