Intraband divergences in third order optical response of 2D systems

The existence of large nonlinear optical coefficients is one of the preconditions for using nonlinear optical materials in nonlinear optical devices. For a crystal, such large coefficients can be achieved by matching photon energies with resonant energies between different bands, and so the details of the crystal band structure play an important role. Here we demonstrate that large third-order nonlinearities can also be generally obtained by a different strategy: As any of the incident frequencies or the sum of any two or three frequencies approaches zero, the doped or excited populations of electronic states lead to divergent contributions in the induced current density. We refer to these as intraband divergences, by analogy with the behavior of Drude conductivity in linear response. Physically, such resonant processes can be associated with a combination of inraband and interband optical transitions. Current-induced second order nonlinearity, coherent current injection, and jerk currents are all related to such divergences, and we find similar divergences in degenerate four wave mixing and cross-phase modulation under certain conditions. These divergences are limited by intraband relaxation parameters, and lead to a large optical response from a high quality sample; we find they are very robust with respect to variations in the details of the band structure. To clearly track all of these effects, we analyze gapped graphene, describing the electrons as massive Dirac fermions; under the relaxation time approximation, we derive analytic expressions for the third order conductivities, and identify the divergences that arise in describing the associated nonlinear phenomena.


I. INTRODUCTION
Motivated by the novel optical properties of graphene, 1-3 many researchers have turned their attention to the linear and nonlinear optical response of 2D systems more generally. 4,5 While there are certainly strong-field excitation circumstances under which a perturbative treatment will fail, 6-9 for many materials a useful first step towards understanding the optical response is the calculation of the conductivities that arise in an expansion of the response of the induced current density in powers of the electric field. 10 In materials where inversion symmetry is present, or its lack can be neglected, the first non-vanishing nonlinear response coefficient in the long-wavelength limit arises at third order, and that is our focus in this paper. The simplest approach one can take to calculate such response coefficients is to treat the electrons in an independent particle approximation, 11 describing any electron-electron scattering effects and interactions with phonons by the introduction of phenomenological relaxation rates. Such a strategy certainly has its limitations, but at least it identifies many of the qualitative features of the optical response, and in particular it identifies what we call "divergences" in that response. We use this term to refer to frequencies or sets of frequencies for which the optical response would diverge in the so-called "clean limit," where the relaxation rates are omitted from the calculation. Under these conditions the predicted magnitude of a response depends critically on the values chosen for the phenomenological relaxation rates. These "divergences" are of particular interest to experimentalists because they indicate situations where the optical response can be expected to be large; they are also of particular interest to theorists since they indicate conditions under which a more sophisticated treatment of scattering within the material, or perhaps a treatment of the response more sophisticated than the perturbative one, is clearly in order.
The optical response of a crystal arises due to interband and intraband transitions. 11 Resonances can be associated with both transitions. For linear optical response, only a single optical transition is involved. A single interband transition can be on resonance for a large range of photon energies, as long as the photon energy is above the band gap. But the resonant electronic states are limited to those with the energy difference matching the photon energy, which depend on the details of the band structure. A single intraband transition can be on resonance only for zero photon energy, and for electronic states at the Fermi surface.
Thus, whether or not these resonances lead to a divergent optical response can differ strongly from material to material. For the nonlinear optical response, intraband and interband transitions can be combined, leading to complicated nonlinear optical transitions. 12,13 As with single intraband transitions, when the nonlinear optical transitions involve the same initial and final electronic states the resonant frequency, which is the sum of all involved frequencies, is also zero. However, due to the interplay of interband transitions, the incident frequencies need not necessarily be all zero for there to be a divergence, and the involved electronic states need not necessarily be around the Fermi surface. Thus, their existence and characteristics are of a more universal nature. It is these divergences that we address in this paper for 2D systems. Although our discussions are illustrated in 2D systems, the underlying physics is the same for systems of different dimension.
Because the nonlinear transitions involve a number of frequencies, these divergences can be classified into different types, associated with different types of nonlinear phenomena.
Several of these types have been widely studied in literature, usually within the context of a particular material or model or excitation condition; the connection with the general nonlinear conductivities is seldom discussed. Our goal here is to demonstrate the universal nature of the expressions for the response, and the generality across a range of materials.
The first type can be called "current-induced second order nonlinearity" (CISNL). It arises when free electrons in the system are driven by an applied DC field; the induced DC current breaks the initial inversion symmetry, and thus generate an effective second-order response to applied optical fields, leading to phenomena such as sum and difference frequency generation. The nature of the divergence here is in the response to the DC field, similar to a single intraband resonance, which would be infinite if phenomenological relaxation terms were not introduced; however, when written as proportional to the induced DC current the effective second-order response coefficients are finite. This phenomenon has been investigated extensively in different materials, both experimentally 14-18 and theoretically [19][20][21] . A second type is "coherent current injection" (CCI), [22][23][24] where the presence of fields at ω and 2ω − δ leads to a divergent DC response as δ → 0 if the excitation at 2ω is able to create free carriers; the divergent response signals the injection of current by the interference of one-photon absorption and degenerate two-photon absorption amplitudes. This is the most widely studied process, both experimentally 25,26 and theoretically. [27][28][29][30] Recent theoretical work has also identified an injection process associated with one-photon absorption and the stimulated Raman process. 28,31 A third type is the jerk current, 32,33 which is a high order divergence involving both a static electric field and an optical field. The static DC field can change the carrier injection rate induced by the optical field, as well as a hydrodynamic acceleration of these optically injected carriers.
We can also identify new divergences, which have not been well recognized in the literature, for two familiar third order nonlinear phenomena. The first arises in cross-phase modulation (XPM) when fields at ω p and ω s are present. The response for the field at ω s due to the field at ω p can diverge when ω p is above or near the energy gap, leading to a phase modulation of the field at ω s that is limited by a relaxation rate. The second also involves excitation with fields at ω p and ω s , but focuses on the degenerate four-wave mixing (DFWM) field generated at 2ω p − ω s . As ω s → ω p this term diverges for ω p above or near the energy gap. These cases merge as ω p → ω s , which corresponds to the most widely studied nonlinear phenomenon of Kerr effects and two-photon absorption. [34][35][36][37][38][39][40][41][42] The very large variation of the extracted values of the nonlinear susceptibilities associated with these phenomena 5,41,43 may be related to such divergences.
In section II we review the general expressions for the third order optical response in the independent particle approximation, and identify in general the divergences that appear associated with the nonlinear optical transitions with a vanishing total frequency. In section III we specialize to the case of gapped graphene, and use it as an example to illustrate the divergences. In section IV we point out the differences between the divergent behavior of gapped and ungapped graphene. In section V we conclude.

II. THE THIRD ORDER RESPONSE COEFFICIENT
The general third order nonlinear susceptibility have been well studied in literature for a cold intrinsic semiconductor, 11 with a large effort devoted to working out many subtle features. In this section we mainly repeat the same procedure for a general band system, and classify the expression in a way that the divergent term can be easily identified.
Writing the electric field E(t) as and other fields similarly, to third order in the electric field the induced current density J (3) (t) is characterized by the response coefficients σ (3);dabc (ω 1 , ω 2 , ω 3 ), where superscripts a, b, ... indicate Cartesian components and are summed over when repeated. The coefficients σ (3);abcd (ω 1 , ω 2 , ω 3 ) can be taken to be symmetric under simultaneous permutation of (bcd) and (ω 1 , ω 2 , ω 3 ), and since the sums over the ω i are over all frequencies there are "degeneracy factors" that arises under certain combinations of frequencies. For example, if fields at ω p and ω s are present we have Often the response coefficient We do not do that here to avoid cluttering the notation, but we consider it implicit in that when we picture these response coefficients we draw arrows associated with all four of the variables appearing in σ (3);dabc (−ω 1 − ω 2 − ω 3 ; ω 1 , ω 2 , ω 3 ), upward arrows associated with positive variables and downward arrows associated with negative variables.
To calculate the response coefficients in the independent particle approximation we label the bands by lower case letters n, m, etc., and the wave vectors in the first Brillouin zone by k. The density operator elements ρ nmk associated with bands n, m and wave vector k satisfy the equation of motion 11 Here we describe the interaction of light with the matter using the "r·E" approach rather than the "p·A" approach involving the vector potential A(t), for the latter can lead to false divergences associated with the violation of sum rules when the number of bands are inevitably truncated to make any calculation. The coefficients ξ nmk are the Berry connections, using the definition in the work by Aversa and Sipe. 11 The interband optical transitions are identified by the off-diagonal terms of ξ nmk , while the rest of terms associated with E(t) are associated with the intraband optical transitions. The last term, ∂ρ nmk ∂t scat , describes the relaxation processes. In our approach, we take a relaxation time approximation, and specify different relaxation time for different transitions, as given below. We solve Eq. (4) perturbatively by setting ρ(t) = j≥0 ρ (j) (t), where ρ (0) (t) = ρ 0 stands for the density matrix in the thermal equilibrium, and ρ (j) (t) ∝ [E] j . The iteration for each order is given by where we take the relaxation terms 44 as nmk , for n = m .
Here the Γ With the evolution of the density matrix, the current density can be obtained through where v mnk are the matrix elements for a velocity operator. A third order response calculation 10 leads to the result with ω = ω 1 + ω 2 + ω 3 and the unsymmetrized coefficients σ (3);dabc (ω, ω 0 , ω 3 ) take the form with Note that the actual value of the energy appearing in, for example, v 0 depends on the corresponding frequency (here ω 0 , a sum of two of the incident frequencies) appearing in The general expression for the conductivity in Eq. (9) immediately indicates the possibilities of the nonlinear phenomena discussed in the Introduction. For CISNL, the conductivities σ (3);dabc (ω 1 , ω 2 , 0) include divergences associated with v 3 → 0; for coherent current injection, the conductivities σ (3);dabc (−2ω, ω, ω) include divergences associated with v → 0; the jerk current is a special case of CISNL, described by σ (3);dabc (ω, −ω, 0), with divergences associated with v 3 → 0 and v → 0; for XPM and DFWM, the conductivities σ (3);dabc (ω, ω p , −ω p ) and σ (3);dabc (−ω s , ω p , ω p ) include divergences associated with v 0 → 0. In special cases, there may be extra divergences identified by a combination of these limits, and the detailed divergences types are determined by the values of S i . Of course, for finite relaxation times we will not have a vanishing v 3 , v 0 , or v. Nonetheless, for frequencies where the real part of one of these quantities vanishes the term(s) in Eq. (9) containing this quantity will make the largest contribution, and we refer to them as the "divergent contributions." Our focus is the identification of these divergent contributions.

III. NONLINEAR OPTICAL CONDUCTIVITY OF A GAPPED GRAPHENE
We apply the approach to a gapped graphene, a two dimensional system. The low energy excitations exist in two valleys, which can be described by a simplified two band model for the unperturbed Hamiltonian The quantity k = k xx + k yŷ is a two dimensional wave vector, ∆ ≥ 0 is a mass parameter to induce a band gap 2∆, τ = ± stands for a valley index, and v F is the Fermi velocity.
Ungapped graphene corresponds to the limit ∆ → 0. All the necessary quantities for the calculation of the third order conductivity are given 45 by For the special two band system we use s = ± to indicate the upper (+) and lower (−) bands.
In this approximation for the dispersion relation of the bands, all the S dabc j can be analytically obtained, and are listed in the Appendix B. Furthermore, the unsymmetrized conductivity σ (3);dabc (ω, ω 0 , ω 3 ) can be written as the sum of two terms, includes all terms in Eq. (9) in which v 3 appears, and σ , and the symmetry {x ↔ y}. We list the independent nonzero components of these tensors, taking S dabc j (· · · ) as an example, as a column vector For gapped graphene, the expressions for S i show features similar to the corresponding expressions for graphene: 44 They are functions of the effective gap parameter E c = max{∆, |µ|} and ∆, in addition to the energies w, w 0 , and w 3 . In all these expressions, E c appears only We can write each S dabc j as a linear combination of these functions, with their arguments themselves being functions of w, w 0 , w 3 , and ∆. All terms can be expanded in even orders of ∆, particularly as functions proportional to ∆ 0 , ∆ 2 , and ∆ 4 , as shown in Appendix B.
With all these analytic expressions in hand, any third order nonlinear conductivity can be calculated and studied directly. In the following, we consider the intraband divergences that are of interest here, by giving the leading contributions to the conductivities.

A. CISNL
We begin by considering the response coefficient σ (3)dabc (ω 1 , ω 2 , 0), where we assume that neither ω 1 nor ω 2 vanishes, and that their sum is finite. This describes a second-order optical nonlinearity induced by a DC field (ω 3 = 0), and if there are free carriers we would expect a divergent contribution because of the large response of the free carriers to the DC field. A picture of the excitation scenario is given in Fig. 1. In our expression for σ (3);dabc (ω, ω 0 , ω 3 ) the divergence of interest is clearly signaled by v 3 → 0, and the divergent response term can be obtained from the analytic expressions. When this is isolated, the most important relaxation parameter is Γ a , which we treat differently than we treat the other relaxation parameters; for j = 2, 3 we set Γ e to give v = w and v 0 = w 0 in the expression for σ (3);dabc . The divergent term can be written as with σ 3 = σ 0 ( v F e) 2 /π, σ 0 = e 2 /(4 ), and We refer to the terms in σ (3)dabc (ω 1 , ω 2 , 0) that are not divergent as v 3 → 0 as the "fieldinduced second order nonlinearity." They exist in the absence of free carriers, and have an analogue in the field induced second order nonlinearity of usual semiconductors, which leads to processes such as electric-field induced second harmonic generation (EFISH). 44,46 Adopting this perspective, we can write the effective second order conductivity σ (2);dab (ω 1 , The first term characterizes the current-induced second order response coefficient, and arises from the divergent contribution in Eq. (18), with all other relaxation times taken as Γ (j) a/e = Γ. All the remaining contributions are collected in the field induced response coefficient σ  In cross-phase modulation the propagation of light at frequency ω s is modified by the presence of an intense optical field at a different frequency, ω p ; it is characterized by the response coefficient σ (3) (ω s , −ω p , ω p ). This conductivity exhibits a divergence as w 0 → 0.
Unlike CISNL, the divergence arises from the nonlinear response of the system and a static field is not required. The process is pictured in Fig. 3; we will see below that the divergent term vanishes unless ω p is near or above the effective band gap 2E c . Cross-phase modulation of a signal frequency ω s can be effected not just by a CW beam at ω p , but by a pulse of light centered at such a frequency. So in general we seek σ (3) (ω s , −ω p + δ, ω p ), where knowledge of this expression for a small range of frequencies ω p centered about a nominal pump frequency, and for a small range of detunings δ centered around zero, will allow us calculated the crossphase modulation of a signal by a pump pulse. As δ, Γ (2) a → 0, the leading term of the relevant unsymmetrized conductivity is (1) , where we set all first order and third order relaxation energies to be the same for simplicity, and and .
Here we used We note that Z 1 and Z 2 themselves diverge as x → −x 3 (i.e., ω p → ω s ), which will be discussed in next section.
Using the expression above we can write where σ xpm;d (ω s , ω p ; δ) contains a divergence with respect to δ + iΓ and σ xpm (ω s , ω p ; δ) contains the non-divergent response. Note that the divergent term is independent of the sequences of the limits δ → 0 and Γ → 0.  and ω p → 0, as would be expected from the discussion of CISNL above, and as ω s → ω p , as we discuss in the section below. As would be expected from the discussion of CISNL above, there are other divergences or resonant peaks associated with either ω s → 0 or interband transitions; the latter are not our focus in this work. Away from these resonances, the values of the conductivity for the cases ω p /E c = 1.5 are much smaller than those for the cases ω p /E c = 2.5, which include the intraband divergences, and the contribution from σ (3);xxxx xpm;d (ω s , ω p ; 0) is generally the largest part of σ (3) (ω s , −ω p , ω p ) away from these other divergences. We can get some insight into the importance of pump frequency by noting that as Γ → 0 the divergent term σ and T (y; x) = G(y; x + i0 + ) + G(y; −x + i0 + ) = 2iπθ(x − 2y).
So in the clean limit σ xpm;d (ω s , ω p ; δ) vanishes if ω p is below the effective gap 2E c . The frequency variation δ comes from the pumping beam, which is negligible for very long pulse duration; then σ (3) xpm;d is a pure imaginary number, the modulation is for the refraction rather than the absorption. Deviations from this arise from the finite relaxation rate for other transition processes, but clearly this divergence is associated with resonant excitation from the valence to conduction band at ω p , as sketched in Fig. 3.

C. DFWM
In Eqs. (23) and (30) there appear to be divergences as v → ±v 3 . Working out the expression in detail it can be immediately seen that in fact no divergence results as v → v 3 , but a divergence does indeed arise as v → −v 3 . The leading divergence is associated with It gives a higher order divergence with a Lorentz type lineshape, as shown in Figs. 4 (b) and (d). Such divergence also exists in a widely studied nonlinear phenomena of four-wave mixing, which is characterized by the response coefficient σ (3) (ω p , ω p , −ω s ). As for XPM, the divergent point is at ω s = ω p , but since in general the generated field is at 2ω p − ω s rather than at ω s the path to the divergence is different.
As ω s → ω p , the divergent term arises in the unsymmetrized conductivities σ (3) (2ω p − ω s , ω p − ω s , −ω s ) and σ (3) (2ω p − ω s , ω p − ω s , ω p ). Taking δ = ω s − ω p and all Γ (j) a/e as small quantities, the conductivity can be approximated as Here we set all the relaxation parameters Γ (j) a/e except Γ (2) a equal to Γ. The exact expression can be obtained following Eq. (22), or even from the full expression in Appendix. B. The physics of this divergence can be revealed if we consider the clean limit for Z 3 and Z 4 . They We illustrate these divergences in Fig. 6 for two different pump frequencies, ω p /E c = 1.5 and 2.5. We plot the results for an undoped system, |µ|/∆ = 0 , and a doped system with µ|/∆ = 1.4. The conductivity strongly depends on the ratio ω/E c . For ω/(2E c ) < 1, the divergent term vanishes in the clean limit and only the non-divergent term survives.
The features in Fig. 6(a) as ω s is around ω p are induced by the nonzero relaxation parameters. For ω/(2E c ) > 1 the divergent term survives and to leading order varies as where we separate out the second order intraband relaxation parameter, and put the rest equal to Γ. The near degenerate FWM is approximately determined by Note that the nonlinear response to a pulse of light can expected to be very complicated due to the strong dependence on the detuning.

D. CCI
Now we turn to another divergence associated with v → 0 in Eq. (9), which exists in the conductivity σ (3) (−2ω + δ, ω, ω) as δ → 0. This divergence describes coherent current injection. Setting all the relaxation parameters Γ leading contribution with respect to the small quantity Γ a and δ given by While the full expressions of η cci (ω, Γ) can be obtained from our general analytic expressions, the underlying physics can be easily shown at the clean limit; setting Γ → 0 + , we find where the functions g 1 and g 2 are given by The coefficient η cci describes the two-color coherent current injection coefficient. In the clean limit, it includes two terms: the one involving g 2 is nonzero only for ω > 2E c , and the one involving g 1 is nonzero for ω > E c . Both terms arises because of the interference of one-photon and two-photon absorption processes. However, the contribution from the term involving g 1 comes from the interference between the pathways of one-photon absorption and degenerate two-photon absorption, as illustrated in the transition process on the left side of Fig. 7; while the term involving g 2 comes from the interference between one-photon absorption and the stimulated Raman process, as shown in the transition process on the right side of Fig. 7. In experiments involving typical semiconductors 2 ω is usually greater than the band gap energy but ω is not; however, if ω also is greater than the gap there can be a contribution due to stimulated electronic Raman scattering. In the clean limit, a direct calculation of the functions g i (x) gives g 1 (1) = g 2 (2) = 0. For the component σ (3);xxxx , their contributions are maximized at x ≈ 1.2 for g 1 (x) and x ≈ 2.4 for g 2 (x). However, at x ≈ 2.4, the contribution from the stimulated electronic Raman scattering is no more than 20% of the contribution from the usual injection process; the total contribution has no maximum around x ≈ 2.4.
A direct consequence of the divergence v → 0 is the form of the current response to a pulse light. For simplicity, we only take the beam at 2ω as a pulse, and then the time evolution of the current density is With substitution of the leading term in Eq. (35) we can get Here E 2ω (t) and E ω (t) are the time evolution of pulses with center frequencies at 2ω and ω, respectively. The first term at the right hand side is a source term, while the second term describes the damping. This equation shows exactly how the injection process occurs.
As an illustration, we plot the σ (3);xxxx (−2ω, ω, ω) as a function of ω for different chemical potentials |µ|/∆ = 0, 1.1, 1.2, and 1.4; the other parameters are taken as Γ a /∆ = 0.01, Γ/∆ = 0.05, and ∆ = 1 eV. In the clean limit, σ (3);xxxx is pure imaginary, and although with the inclusion of damping the real parts do not vanish, they are about an order of magnitude smaller than the imaginary parts, as shown in Fig. 8. In the calculations at finite relaxation parameters, both the real and the imaginary parts show obvious peaks/valleys around ω ∼ E c and ω ∼ 2E c , which correspond to the contributions from the terms including g 1 and g 2 , respectively. In contrast to the situation discussed after Eq. (38) for clean limits, the appearance of the peaks for ω > 2E c arises because of the inclusion of the finite relaxation parameters. There also exists increases of the conductivity values as ω → 0, mostly for nonzero µ/∆. They are associated with a divergence induced by the free-carriers, which is not our focus in this work.
In general, the leading orders are The clean limits of Q 1 and Q 2 are with α = E c /∆. Here Q 1 (x; 0) is nonzero only when the one-photon absorption exists as ω > 2E c , while Q 2 (x; 0) exists only for a doped system where α > 1. These two terms have a different power dependence, and the term involving Q 2 (x; 0) can exists for any optical field frequency. Thus the frequency dependence of the response can be used to distinguish between them.

IV. COMPARISON WITH GRAPHENE
Some of the phenomena discussed above have been considered earlier for doped graphene. 19,44,47 Although the expressions we derived above are normalized to the gap parameter ∆, it is safe to take the limit as ∆ → 0. This is because our general expressions of the conductivity can be safely reduced to the case of graphene with taking ∆ = 0, as shown in Appendix B.
The results of CISNL, DFWM, and CCI in graphene have been discussed earlier 19,44,47 in the clean limit and with finite relaxation parameters. Here we give a brief discussion for XPM and the jerk current in graphene. In the clean limit, the XPM for graphene can be found from Eq. (29) to be Here the chemical potential induced gap plays a role similar to that of the gap parameter in gapped graphene. For the jerk current, the leading term becomes

V. CONCLUSIONS
We have systematically discussed intraband divergences in the third order optical response, and identified the leading terms in the corresponding third order conductivity. Due ot the combination of intraband and interband transitions, these divergences can appear at optical frequencies, and lead to large nonlinear conductivities. We have shown that the existence of such divergences is very general, independent of the details of the band structure.
We illustrated these divergences in gapped graphene, with third order conductivity obtained in an analytic expression for massive Dirac electrons under relaxation time approximation.
Such divergences are of interest to experimentalists, because within the independent particle treatment presented here the optical response is limited only by the phenomenological relaxation times introduced in the theory, and thus that optical response can be expected to be large. As well, at a qualitative level the predicted nature of the divergent behavior is robust against approximations made in describing the details of the interband transitions.
The divergences are also of interest to theorists, because one can expect that under such conditions the kind of treatment presented here is too naive. This could be both because more realistic treatments of relaxation processes are required, and as well because the large optical response predicted could be an indication that in a real experimental scenario the perturbative approach itself is insufficient. Thus the identification of these divergences iden- with ω = ω 1 + ω 2 + ω 3 and ω 0 = ω 2 + ω 3 . Using the iteration in Eq.
It shows that the diagonal terms of the Berry connections ξ nnk appear always with the derivative ∇ k to form a gauge invariant term. 11 The third order terms are