Anomalous optical saturation of low-energy Dirac states in graphene and its implication for nonlinear optics

We reveal that optical saturation of the low-energy states takes place in graphene for arbitrarily weak electromagnetic fields. This effect originates from the diverging field-induced interband coupling at the Dirac point. Using semiconductor Bloch equations to model the electronic dynamics of graphene, we argue that the charge carriers undergo ultrafast Rabi oscillations leading to the anomalous saturation effect. The theory is complemented by a many-body study of the carrier relaxations dynamics in graphene. It will be demonstrated that the carrier relaxation dynamics is slow around the Dirac point, which in turn leads to a more pronounced saturation. The implications of this effect for the nonlinear optics of graphene are then discussed. Our analysis shows that the conventional perturbative treatment of the nonlinear optics, i.e. expanding the polarization field in a Taylor series of the electric field, is problematic for graphene, in particular at small Fermi levels and large field amplitudes.

Graphene is a two-dimensional material made of carbon atoms in a honeycomb structure.Its reduced dimensionality and the symmetries of its crystalline structure render graphene a gapless semiconductor [1].Graphene exhibits a wealth of exceptional properties, including a remarkably high mobility at room temperature [2], Klein tunneling and Zitterbewegung [3,4], existence of a nonzero Berry phase, anomalous quantum Hall effect [5][6][7], quantum-limited intrinsic conductivity [8], and a unique Landau level structure [9,10].Underlying these peculiar electronic properties are its pseudo-relativistic quasiparticles that obey the massless Dirac equation [1].As a direct consequence of their massless nature, the Dirac fermions have definite chiralities [11,12].Owing to the specific symmetries of the crystalline structure of graphene, the dynamics of the massless Dirac quasiparticles and their chiral character are topologically preserved-i.e., many-body induced band renormalizations as well as any moderate perturbations of the lattice will not open a gap in graphene's band structure [13].A large number of the unusual properties of graphene are associated with the topologically protected band-crossing and the chiral dynamics of the charge carriers [3].
One major consequence of the topologically protected chirality of the charge carriers is the anomalous structure of the interband coupling mediated by an electromagnetic field.Its dipole matrix element obtained in the length gauge [14] exhibits a singularity at the degeneracy points, in contrast to ordinary (and even other gapless) semiconductors [15,16].This has raised some controversy regarding the treatment of the optical response of graphene [15,17,18].Specifically, the perturbative treatment of the nonlinear optical response has been questioned [17,19].The nonlinear optical coeffi-cients of graphene obtained by means of perturbation theory suffer from a nonresolvable singularity [15,17].Although substantial effort has been spent on developing comprehensive models for the nonlinear optical response of graphene [17,[19][20][21][22][23][24], a self-consistent theoretical model that can resolve the above issue is still lacking.In addition, many experimental studies of the nonlinear optics of graphene have been reported [25][26][27][28][29][30]-some of these studies are difficult to reconcile with existing theoretical models.
In this Letter, we show that the singular nature of the interband dipole coupling has some significant physical implications: it causes the charge carriers in the vicinity of the Dirac points to undergo ultrafast Rabi oscillations accompanied by slow relaxation dynamics, which, intriguingly, yields an anomalous saturation effect.This finding necessitates revisiting the perturbative treatment of the nonlinear optical response of graphene to account for the extreme nonlinear interactions around the Dirac points.These conclusions will be reached by describing the dynamics of the charge carriers with semiconductor Bloch equations (SBEs) [31,32].
We consider a free-standing graphene monolayer (in the xy plane) illuminated by a normally incident electromagnetic field.The monochromatic and spatially uniform optical field at the graphene layer is described by E(t) = E 0 e iωt + c.c., where E 0 is parallel to the graphene plane.The light-matter interaction is considered semiclassically and the external field coupling is obtained in the length gauge [14].For photon energies below approximately 2eV, the electronic dynamics of the quasiparticles in the absence of external radiation is adequately described by the massless Dirac equation yielding the relativistic energy-momentum dispersion E k = ± v F |k| [1], where k is the Bloch wave vector with respect to the Dirac point and v F is the Fermi velocity.
The SBEs describe the coupled dynamics of the population difference N (k, t) and the polarization (coherence) P(k, t) in the momentum state k.In the absence of electromagnetic radiation, the population difference relaxes to , where f (E) is the Fermi-Dirac distribution.An electromagnetic field drives the system out of equilibrium via the coupled intraband and interband dynamics.In a moving frame {τ, k } = {t, k − δk(t)}, where δk obeys ∂δk ∂t + Γδk = − e E(t) (Γ is a phenomenological intraband relaxation coefficient), the dynamics of the charge carriers is governed by where Φ(k, t) = eE• φk k is the matrix element of the external potential of the direct optical transition, and the unit vector φk is defined as φk = ẑ × k/k.The frequency k = 2E k is the energy difference between the energy levels of the conduction and valence bands.γ k and γ (2) k are k-dependent relaxation coefficients stemming from many-body effects such as electron-electron and electronphonon interactions.It is worth pointing out that the moving frame accounts for the coherent Bloch oscillations, which play an important role in high-harmonics generation [33].The detailed derivation of our theoretical model is provided in the Supplemental Material [34].
The light-graphene interaction as described by Eqs.(19a)-(19b) can be interpreted as an ensemble of inhomogeneously broadened two-level systems (one for each k).The last term in each of the two equations will lead to Rabi oscillations.Because of the singularity in Φ(k , τ ) for |k| → 0, we can expect ultrafast Rabi oscillations around the Dirac point, which are damped by many-body interactions.The decay terms drive the twolevel systems towards an equilibrium state.Since the interband coupling is strong around the Dirac point (equivalent to highly intense illumination), the effective field leaves the two-level systems in a statistical mixture of the ground and excited states with equal weights and absorption quenching takes place.Thus, the states around the Dirac points undergo a saturation effect, even when illuminated by an arbitrarily weak electromagnetic field.
This saturation behavior can be further understood by studying the steady-state solution of the SBEs, which is where Φ k = eE 0 • φk / k is the complex phasor associated with Φ(k, t).The function ∆ k = ω − k denotes the detuning of the two-level system at k with respect to the excitation.Since | Φ k | is arbitrarily large for small-k states, in the vicinity of the Dirac point, the population N k cannot be expanded in a Taylor series of the field Φ.This implies that the nonlinear optics of graphene is in principle a nonperturbative problem.Indeed, due to the singularity of the interband coupling in graphene, there is always a region around the Dirac point where graphene is optically saturated.The saturation threshold E sat k is given by k . ( Saturation occurs, of course, in any two-level system at high field intensities.However, in graphene the saturation threshold field E sat k is zero at the Dirac point (k → 0) and, hence, there is always a region of k-space where E 0 > E sat k , even for arbitrarily weak intensities.The peculiar low-threshold saturation mechanism in graphene can be quantitatively resolved using a timedomain analysis of the graphene SBEs.For the sake of comparison, this analysis has been performed for two distinct continuous excitations with optical frequency of ω = 80meV (terahertz range) and ω = 800meV (infrared), respectively.In both cases, the electric field is linearly polarized along the ŷ direction with magnitude E 0 = 10 6 V/m.Graphene is assumed to be undoped here and is initially held at room temperature.The relaxation coefficients γ (1) k and γ (2) k are determined using a microscopic theory, which encompasses carriercarrier as well as carrier-phonon scattering channels and takes into account all relevant relaxation paths including interband and intraband and even inter-valley processes [36,37].The reader is referred to the Supplemental Material for the details of the many-body model [34] as well as the methodology used to extract the relaxation coefficients.The resulting relaxation coefficients are plotted in Fig. 1(b) and (c).We note in particular that γ The relative change in the stationary component of the population difference due to the optical excitation as well as the amplitude of the oscillating induced polarization are shown in Fig. 1(d) and (e).To obtain the steadystate components, we performed Fourier analysis within a time window where the transient response has died out.As expected, a well-pronounced modified population difference around the Dirac point due to the spontaneous polarization effect (dark red region around the center) is observed.This effect is stronger for lower-frequency excitations-indeed, according to Eq. (3) a smaller detuning yields a weaker saturation threshold.The region in k-space where the spontaneous optical saturation is significant is well extended from the Dirac point.We note here that the size of the region depends on the applied field intensity-we will show below that this is the origin of the nonperturbative nature of the nonlinear optical response.In addition, there is, of course, the traditional optical saturation region for ∆ k ≈ 0 (indicated by the yellow dashed line).
Before we continue with the importance of this anomalous optical saturation for the nonlinear optics of graphene, let us briefly make a few remarks.First, the origin of the inverse dependence of the interband transition matrix element on the wavenumber can be linked to the distinctive mathematical structure of the current operator.It is straightforward to show that the interband coupling matrix element at wavenumber k is rcv where J cv is the offdiagonal element of the current operator.In contrast to ordinary semiconductors, the off-diagonal components of the current operator in graphene and other chiral materials are strictly nonzero even at the band crossing points [38].As a direct consequence of this property of massless Dirac quasiparticles, the interband part of the position operator carries a first-order singularity at the degeneracy point.Second, one may wonder why we have not used the velocity gauge, in which optical coupling is obtained by minimal substitution k → k + eA where E = −∂A/∂t.It is important to note that this approach is not gauge invariant in the "effective Hamiltonian" picture [39][40][41].We show in the Supplemental Material that a modification of the velocity gauge is indeed required to yield a physically correct result [34].This modification gives rise to the 1/k dependence of the interband coupling in the vicinity of the Dirac point [41].Third, although SBEs can model the essential physics, a number of approximations have evidently been made in deriving them.For instance, for intense sub-terahertz excitations, quasi-instantaneous thermal effects can cause a population pulsation at a rate faster than the interband Rabi oscillations and thus they cannot be modelled by a spectral broadening [35].However, owing to the fact that the anomalous saturation effect occurs around the Dirac point where the quasiparticles are off-tuned with the excitation photons, quasi-instantaneous thermal effects-taking place dominantly around the zerodetuning circle-have minimal impact on the anomalous saturation effect.Moreover, since around the Dirac point population difference is fairly independent of temperate T , i.e., ∂N eq k /∂T ≈ 0, temporal fluctuations of the temperature will not significantly obscure the region where saturation happens.
We now turn to the nonlinear optics of graphene.Let us consider a nonlinear pump-probe experiment in which graphene is simultaneously subjected to the a pump (ω c ) and a weaker probe (ω p ) laser beam.The conductivity tensor of graphene in the presence of the pump field and "seen" by the probe field is calculated in the Supplementary Material [34].Fig. 10 displays the change to the conductivity tensor due to the pump fields described earlier.The real part of the conductivity is re- lated to the absorption coefficient of the probe beam via α ≈ Re{σ}/4ε 0 c, where c is the speed of light in vacuum.The relative change in absorption of the probe beam is also shown in Fig. 10.For the 80meV pump beam, there is strong saturation for a probe beam at the same frequency-this is because the pump beam has saturated the interband transition.For a probe beam with a much lower frequency (ω p ≈ ω c /100), there is also strong saturation-which must be due to the unconventional effect discussed above.For the 800meV pump, a weaker saturation is observed for a low-frequency probe.Indeed, for a higher-frequency pump the region of the low-threshold saturation is smaller [observe in Eq. ( 3) that a larger detuning results in a larger saturation field].Although the above example discusses the anomalous saturation effect for undoped graphene, the effect happens for a range of Fermi levels.Therefore, it will also be observed in samples where the charge-neutrality point fluctuates in space.
We now move on to discuss the applicability of perturbation theory in the analysis of the nonlinear response of graphene.As detailed in Ref. [15], the standard perturbative treatment of the optical response of graphene leads to a nonresolvable singularity in its higher-order nonlinear coefficients (beyond the linear response) originating from the small-k states.We now understand that these states should be excluded because they are saturated.However, since the saturation threshold depends on the field intensity, the nonlinear response coefficients calculated with standard perturbation theory become field dependent too.But let us investigate under what circumstances the standard nonlinear coefficients have meaning.
As mentioned earlier, neglecting the momentum of the absorbed photon, the optical transitions are vertical in k-space and, therefore, every point in the reciprocal space can be treated independently.It is argued in Refs.[15,17,22] that nonlinear frequency mixing in graphene can be decomposed into a number of additive contributions (see the Supplemental Material for a complete theoretical analysis [34]), i.e., the nonlinear conductivity tensor is σ (l) ≈ k,|k|>Ksat I (l) k , where K sat is the radius (with respect to the Dirac point) of the spontaneously saturated region in k-space and is obtained from Eq. ( 3): γ (1) , (4) and k represents the contribution of the quasiparticles with Bloch index k to the l'th order nonlinear optical response.E 0 is the magnitude of the largest electric field component (most often a pump field) participating in the nonlinear process and ω is its frequency.
In order to gain insight into the intensity dependence of the nonlinear response coefficients obtained from the above-described semiperturbative approach, we compare in Fig. 3 the third-order nonlinear response defined as |σ (3) xxxx (ω, ω, −ω)| (where we have now used kindependent relaxation constants as usual) with the results of the full solution of SBEs (see the Supplemental Material [34]).The yellow shaded curves display the nonperturbative solution and the black dotted curves are the results of semiperturbative approach.First, owing to the low-threshold saturation effect, there is a noticeable field dependence of the third-order nonlinear (Kerr-like) response for lower Fermi energies.When the field intensity becomes large enough to extend the saturation region to the excited k-states, the semiperturbative approach fails.As the Fermi energy becomes larger, the optically induced Pauli blocking becomes less important as the low-energy states are already Pauli blocked.It is worth pointing out that we have observed significant dependence of the results on K sat .Therefore, the exact exclusion of the saturated region is necessary to achieve accurate results.
In conclusion, we have demonstrated that the topologically protected singular interband coupling in graphene leads to ultrafast Rabi oscillations, exciting the quasiparticles faster than they can relax back to the ground state.This leads to an anomalous optical saturation of the low-energy quasiparticles in graphene.Subsequently, we have shown that due to this effect the small-k states have to be excluded for the perturbative calculation of the nonlinear optical coefficients of graphene.As a result, the nonlinear coefficients obtained from perturbation theory exhibit noticeable field dependence, particularly for small Fermi levels.Although the effect revealed in this Letter has not yet been observed directly, several experiments have demonstrated the nonperturbative nature of the nonlinear optical response in graphene.For instance, Refs.[25] and [30], which present experimental results for the Kerr nonlinear coefficient of graphene obtained from z-scan measurements, demonstrate that the effective Kerr coefficient is not independent of the intensity of light.At very intense illuminations, the field dependence of the Kerr coefficient might of course have multiple origins, but this effect has been observed even in the weak-field regime, e.g., in Ref. [33,48], which report recent experimental observations of optical and terahertz high-harmonic generation in doped and nearly-undoped graphene.Our theory provides an explanation for these experimental results indicating the nonperturbative nature of the nonlinear optics of graphene.We speculate that similar effects may be found in other Dirac materials and in Weyl semimetals.

Equations of Motion
To begin, we use the low-energy effective Hamiltonian of graphene, which is described by H 0 = v F k • σ in the sublattice (pseudospin) basis.Pauli matrices σ = σi xi are used to expand the operators.For the time being, we ignore the many-body effects such as Coulomb interactions and their collective influence will be phenomenologically included.The calculations are carried out within the length gauge, using the additive potential V ex = eE • r to couple the electromagnetic field to the dynamical equations.The time evolution of the density matrix is then obtained through the Liouville equation as [15] i The 2 × 2 pseudospin density matrix ρk can be expanded in terms of Pauli matrices Substituting Eq. ( 6) into Eq.( 5), one obtains decoupled equations for coefficient n k and vector m k (see Ref. [42]): We will also need the current operator to construct the optical response functions and the current density becomes Having derived the equations of motion in the sublattice basis, we can now switch to the energy diagonal basis.We use "∼" to denote the matrix representation of the operators in the valence and conduction basis.In the energy diagonal basis, matrices 1 0 T and 0 1 T are adopted for the upper and the lower energy levels, respectively.In the energy diagonal basis, the density matrix and the current operator become where k and φk are the unit vectors shown in Fig. 4. At thermal equilibrium, before switching on the incident field, the density distribution obeys Fermi-Dirac statistics where ξkv and ξkc are the many-body annihilation operators for the valence and conduction bands, respectively.The subscript 0 denotes the equilibrium state and E(k) = v F k is the upper energy level.The distribution f (E) is the Fermi-Dirac distribution function where T and µ are, respectively, the temperature and the chemical potential associated with the Fermi energy level E f .

Semiconductor Bloch Equations
It can be shown that the optical response under continuous excitation wave can be estimated based solely on the dynamics of (i) the microscopic population difference N (k, t) and (ii) the microscopic polarization P(k, t) defined as Taking N and P as dynamical variables and using Eq. ( 8), we obtain the equations of motion for the population difference and the microscopic polarization.where Φ(k, t) is essentially the matrix element of the external potential for the given Bloch momentum k describing the direct optical transition between the upper and lower energy levels.This is due to the fact that the momentum of the light is assumed to be negligibly small.
The frequency k = 2E k is the energy difference between the upper and lower energy levels.The coupled equations given in Eq. ( 16) are called the semiconductor Bloch equations (SBEs) for graphene.These equations must be solved simultaneously, and subject to the initial condition imposed by the Fermi-Dirac distribution before turning on the field:

PERTURBATIVE SOLUTION OF THE SEMICONDUCTOR BLOCH EQUATIONS
The SBEs describe the quasiclassical transport and interband excitation problems simultaneously.To convert these dynamical equations into a more convenient form, the transport and interband evolutions are decoupled.The cooperative intraband and interband dynamics can be split by introducing a moving frame in the reciprocal space whose movement is governed by the Boltzmann transport equation [15].The purely classical part of the dynamics drives the distribution function along a trajectory determined by the direction of the electric field.The time-momentum coordinate in the moving frame is designated by {τ, k }.The primed frame is then related to the original frame by {τ, k } = {t, k − δk(t)} where δk obeys where the extrinsic fitting factor Γ accounts for impurity scattering as well as radiation.In the primed frame, the system evolves by pure interband dynamics as described by equations resembling the optical Bloch equations for a generic two-level problem [43]: In compliance with the standard notation of Bloch equations, we introduce w k , u k , v k as where w k is the population inversion in the moving frame.The functions ξ k Φ(k , τ ) and ω k k are also defined for mathematical convenience and are the analytical functions of the exciting field.The function ξ k is the equivalent dipole in the moving frame.The optical Bloch equations describe the interband dynamics for the ensemble of the two-level Bloch states.However, due to the motion of the frame, the intraband dynamics are intertwined with the interband dynamics.The combined nature of the dynamics manifests itself in the time dependence of index k .Since the time scale associated with intraband dynamics is quite different from that of the interband evolution, an adiabatic argument could be applied; however, although the interband transitions acquire an extra time dependence due to the moving frame, they nevertheless take place independently.The coupled Bloch equations can be converted to the standard optical Bloch equations using the two-level approximation [43].Henceforth, we drop the prime in the uvw-coordinate system: where 'dot' designates the time derivative.The function w eq k is the population difference at equilibrium, i.e., w eq k = N eq k .For a weak pump field, the inversion w k tends to relax to w eq k .The coherent terms are the oscillatory functions of the field.To proceed further, the current response in the reciprocal space must be identified.According to Eq. ( 11), together with Eq. ( 12) the induced current is Therefore, the equation of motion describing the time evolution of v k provides enough information to model the interband response of graphene.Neglecting the time variation of ω k in the adiabatic approximation yields Since ω 2 k is much larger than γ 2 2 , we can drop γ 2 2 v k to obtain the result: This set of equations can be solved iteratively, leading to a series expansion for v k and w k .The mathematical structure of the coupled damped-driven harmonic oscillator [Eqs.(24) and (21a)] implies that v k contains only odd powers of the field ξ k and w k contains only even powers of the field: The n'th-order expansion terms w The iterative procedure can be conveniently carried out by defining the operators V k and W k .
In addition to the purely interband multi-photon process described above, a part of the nonlinearity originates from the quasiclassical transport or intraband transitions.As will be clarified further in subsequent sections, the frequency mixing effects in graphene arise from the pure intraband, pure interband and interband-intraband transitions.The pure intraband response occurs only due to the change in the population difference.Using Eq. ( 22), the intraband contribution to the current is then, the n'th-order nonlinearity due to the pure intraband process can be obtained using the following Taylor expansion: For the sake of mathematical convenience, the derivative operator is represented as where the Drude-like coefficient 1/(iω + Γ) is obtained from Eq. ( 18).An intuitive symmetry argument shows that in graphene-a centrosymmetric crystal-even orders of nonlinearity do not exist.The derivations of the linear and third-order conductivity tensors of graphene are discussed below.
Linear Optical Response of Graphene Equation ( 27) for n = 1 gives the intraband conductivity tensor in the k-space where we have assumed that the electric field is E(t) = E p exp(iω p t) + c.c and σ intra (k, ω p ) is defined via intra (k, ω p ) • E Performing the integration in the reciprocal space, the off-diagonal terms vanish and we arrive at where g s and g v are the spin and valley degeneracy factors, respectively.The interband linear optical response of graphene can be obtained using the master equation (24).The linear optical response is a single-photon process and unlike higher-order terms it can be obtained independently for interand intraband contributions.The first-order solution of Eq. ( 24) can be derived by replacing w k and ẇk with w eq k and 0, respectively, Integration over the reciprocal space and including the density of states gives Third Order Frequency Mixing in Graphene Throughout this section we assume that three complex fields with a time dependence of e iωpt , e iωqt and e iωrt are mixing through the third-order conductivity of graphene.As mentioned earlier, the third-order optical response can be interpreted as a three-photon process and different terms contribute to the third-order conductivity tensor, namely a pure intraband term, a pure interband term, and a combination of both.The distinct photon processes contributing to the third-order optics of graphene are schematically shown in Fig. 5.The intraband dynamics causes the quasiparticles to travel along the trajectory determined by the direction of the electric field at the graphene plane.The quasiclassical Boltzman-like dynamics are pictured schematically by displacement of the entire distribution of the quasiparticles over the reciprocal space.The interband contributions are shown by two-level transitions of the quasiparticles predominantly around the zero detuning region.The adopted mathematical structure allows us to easily find the conductivity tensors associated with the six processes shown in Fig. 5.

Pure intraband third-order nonlinearity
The pure intraband contribution is fundamentally part of the third-order nonlinearity which occurs exclusively due to the Boltzman-like transport.The third-order intraband process is schematically displayed in Fig. 5a.Equation (28) determines the nonlinear contribution of the third-order intraband evolution as k,1 (ω p , ω q , ω r ) = −ev F P I k Dk (ω r ) Dk (ω q ) Dk (ω p )N eq k (34) where the normalized gradient operator Dk was introduced in Eq. ( 29).The subscript '1' is used to identify the intraband contribution and the superscript '(3)' refers to the third-order effect.Here we have made use of the intrinsic permutation operator P I ; all possible permutations of the input frequencies ω p , ω q , and ω r must be summed over.The overall intraband conductivity tensor is then obtained as Pure interband third-order nonlinearity Pure interband third-order nonlinearity can be obtained using the master equations ( 24) and (21a) in the moving frame.Power expansion of the inversion and the polarization in terms of the exciting field as in Eqs.(25a) and (25b) gives v  k can be found from the following equations These equations can be solved simultaneously to find the contribution of the interband dynamics associated with Bloch index k to the third-order optical nonlinearity.Making use of the operators V k and W k allows a compact solution as Performing the integration in the reciprocal space and applying the permutation operator yields the final expression for the interband conductivity tensor.The interband conductivity tensor is then obtained as where the summation goes over all quantum states.This term is obviously singular at the origin.

Interband-intraband third-order nonlinearity
The master equations together with the notion of a moving frame lead to four possible combinations of the interand intraband dynamics shown in Fig. 5c-f.
• Interband-Interband-Intraband: According to Fig. 5c, the interband dynamics modify the population difference through a non-parametric transition.The third photon (which is actually the probing photon!) causes an intraband transition.The master equations ( 24) and (21a) yield the modified population accomplished via two subsequent photon processes.
k (ω p ) The third process causes the modified population to move over the reciprocal space and creates an intraband current described by equation (27).The tensor associated with this process reads • Intraband-Intraband-Interband: Following two subsequent intraband transitions at the frequencies ω p and ω q , the population difference oscillates at the sum frequency ω p +ω q .The third photon probes the graphene whose quasiparticles are oscillating over the reciprocal space and the current is induced due to interband transition.
This process is sketched in Fig. 5d and is mathematically encoded by I k,4 (ω p , ω q , ω r ) = −ev F P I φk V k (ω p + ω q + ω r ) Dk (ω q ) Dk (ω p )N eq k (41) • Intraband-Interband-Intraband: Fig. 5e displays the ordering of the processes.Following an intraband transition, the population difference oscillates at frequency ω p .The second process is interband which creates a coherence (polarization) at the frequency of ω p + ω q .The induced polarization is driven by an additional intraband transition to create a current at the frequency of ω p + ω q + ω r .Although the last transition is intraband, the induced current is of interband nature and is modulated by the moving frame.
• Interband-Intraband-Intraband: An interband transition caused by the first photon creates the polarization oscillating at frequency ω p .Due to the Boltzman-like transport, the induced polarization is modulated by two consecutive harmonics ω q and ω r .The induced current has an interband nature and is shown in Fig. 5f.The overall dynamics are encoded by I k,6 : The overall intraband-interband conductivity tensor is then obtained as σ intra−inter (ω p , ω q , ω r ) =

LENGTH GAUGE VERSUS VELOCITY GAUGE IN GRAPHENE
In the long-wavelength regime where the spatial dependence of the electromagnetic field can be neglected, there are two theoretical approaches to couple electromagnetic radiation into the Hamiltonian.Within the so-called velocity gauge, light and matter are coupled via the minimal substitution of the vector potential A(t), with the electric field given by E(t) = −∂A/∂t.The alternative approach is referred to as the length gauge [44] where the field is directly coupled by means of the additive scalar potential V (r) = eE • r with the elementary charge e > 0. Although the velocity gauge preserves translation symmetry and different Bloch states remain uncoupled, there are several undesirable features associated with the velocity gauge that plague the calculations [39,44].
In particular, the treatment of the nonlinear optical response within the velocity gauge is susceptible to numerical errors caused by truncation of the band space [41, 44? ].Since the electric field is proportional to the time derivative of the vector potential, poles at ω = 0 inevitably appear.Diverging terms are not obviously real in semiconductors.It has been rigorously proven that they essentially vanish due to time reversal symmetry and an effective mass sum rule [44].However, a full calculation of the optical transitions over the entire band is required to eliminate the diverging terms [39] and numerical errors due to truncation amplify.Moreover, in two-level systems where an effective Hamiltonian is intended to describe the electrons' dynamics (as in the case of graphene), the application of the velocity gauge is quite vulnerable [41].References [41] and [44] show that in the perturbative treatment of the nonlinear optical response of such systems, the contribution made by the remaining parts of the band cannot be ignored in the calculations.This observation implies that a local interpretation of the optical transitions within the velocity gauge is not reliable and that all Bloch states collectively impact the optical response [41].
To cope with this situation, the length gauge can be employed.The different contributions of the position operators have been extensively discussed in the literature [39,44] and it has been shown that the perturbative treatment of the nonlinear optical response of semiconductors can be reliably performed in a two-level model by making use of the length gauge [45].This gauge eliminates the nonphysical diverging terms and renders it possible to interpret the optical transition locally over the band space.The matrix element of the position operator between the Bloch states indexed by (k, s) and (k , s ) is [39] The position operator can be interpreted as the generator of translation in the space of the Bloch functions [39].The Berry potential ζ ss (k) is required as the geometric phase correction and, therefore, its action in nontrivial topologies is crucial.
To shed light on the applicability of the velocity gauge to perturbation theory, as a benchmark, we proceed with finding the linear optical response of graphene within the velocity gauge.Assume that a graphene layer lying in the xy plane is illuminated by an obliquely incident electromagnetic field E(r) = E 0 exp (iωt − ik 0 .r),where k 0 is the wavenumber of the incident field.For mathematical convenience, the tangential component of k 0 is denoted by q k 0 • (xx + ŷ ŷ).Within the velocity gauge, the time derivative of the divergence-less vector potential A is linearly related to the electric field.The vector potential is The interaction Hamiltonian is then ĤI = A • ˆ J where ˆ J is the current operator derived in Eq. (12).Using the customary minimal substitution prescription, the overal Hamiltonian reads For compactness, we have defined the operator The linear variation of the density matrix due to the presence of the external potential is calculated using the perturbation expansion The first-order perturbation theory gives δ ρ(1) , which linearly depends on the electric field: where f (E) is the Fermi-Dirac distribution.The index s refers to the upper and lower energy states.The parameter γ appearing in the denominator is the phenomenological relaxation coefficient.The first-order induced current is then obtained as J q = Tr δ ρ(1) ˆ R(q) (51) From now on, we define Π ss kk as The unit vectors û and φu are defined as where k (k, q) = (k + q)/ |k + q|.The integrand in Eq. ( 52) is explicitly written as ξ k,q defined as Using Eq. ( 46), the conductivity tensor associated with the transitions between the quasiparticles with Bloch indices k and k+q is then obtained as σ k,q = iξ k,q (ω)/ω.However, as we expected, this term is diverging at lower frequencies due to the pole at ω = 0.It is at this point that the undesirable features of the velocity gauge manifest themselves.The diverging term appears in the optical response and, therefore, special care must be taken.As discussed, this singularity would be resolved in a full-band calculation where all contributions are included.In order to eliminate this singularity, one can manually add a zero to the optical response to cancel out the pole at ω = 0 and yield a physically correct result: Accordingly, we define Πss kk as The conductivity tensor now reads where E i = v F k i is defined to make the integral dimensionless.In the long-wavelength limit, the optical conductivity obtained above and the results of the calculations within the length gauge offer identical expressions for the linear conductivity tensor.It should be noted that where γ ss = Γ is, by definition, the intraband relaxation coefficient.This part of the conductivity is obviously responsible for the intraband dynamics manifested in the ∂/∂k terms appearing in Eq. (45).Likewise, for the interband contribution where s = s where E k,s = −E k,s .The prefactor 1/2E s,k together with the multiplicative vector φk corresponds to ζ ss (k) appearing in Eq. ( 45).
In conclusion, by removing the artificial diverging pole in the velocity gauge, both approaches yield identical results.The aforementioned pole arises due to the band space truncation and it can be removed by developing the sum rule reported in Ref. 44.In the effective two-band model, the velocity gauge should be repaired to account for the nonphysical terms and as a result, the perturbative treatment of the optical response of graphene within the velocity gauge, in its original form, is not perfectly reliable.The length gauge can be consistently employed to develop a full theoretical model for the nonlinear optics of graphene.(II) Use a low-pass filter to get rid of higher harmonics in the polarization envelope as well as the population pulsations.
(III) Curve-fit the decaying part of the dynamics to an exponential function.The decay of the microscopic polarization is γ k and that of the population is γ k .
Figures 6 and 7 display the steps followed to extract the microscopic relaxation coefficients for the pulse excitation of ω = 80meV and 800meV, respectively.It is assumed that the amplitude of the electric fields for the both optical pulses are E 0 = 10 6 V/m.The duration of the pulses have been selected to be 800fs and 400fs for the cases of 80meV and 800meV, respectively.The results of the calculations for a few distinct points over the k-space are shown in Fig. 8.As expected, the low-energy excitation offers significantly slowed relaxation dynamics.It is also noting that, unlike the coherence dephasing coefficient γ    where g s and g v are the spin and valley degeneracy factors, respectively, and D = 1/(2π) 2 is the density of the states in the reciprocal space.The integral goes over the reduced reciprocal space where the Dirac dispersion is valid.A more rigorous treatment would require including additional contributions accounting for the nonlinear frequency mixing due to the pure intraband process.However, the above formulation is adequately accurate for large enough pump frequencies.The interband conductivity consists of two contributions, namely incoherent and coherent terms [47].
Similarly to the intraband part, the incoherent contribution simply results from reduction of the population difference over the reciprocal space.The coherent term, however, enters due to population pulsation at the beat frequency ω c − ω p .It is argued that the coherent term involves interference between the pump and probe fields [47].The population pulsation followed by absorption of a second photon from the pump field acts mathematically as an additional complex-valued contribution to the population difference denoted by δN puls .The overall pump induced interband conductivity tensor is The validity of the theory requires that the coherence length of the pump laser be much larger than the wavelength of the probe field.

NONPERTURBATIVE KERR-TYPE NONLINEARITY OF GRAPHENE
A nonperturbative formulation of the Kerr effect in graphene is obtained by taking the steady-state population difference as the optically modified inversion.In the steady-state analysis, essentially only one-photon processes are retained.The impact of higher-order effects including two-photon absorption (TPA) is then incorporated into the response function via additional complex contributions to the population difference.By substituting the effective complex population difference into the linear response theory, one obtains the induced interband nonlinear current as where the complex function δN T P A k accounts for the population oscillations at the second harmonic of the excitation.In the dressed-state picture, TPA is conceived as the absorption of two subsequent photons via a virtual state.In the density-matrix formulation, that contribution results in the oscillation of the population difference at the second harmonic.By expanding the inversion as N k ≈ N st k + (N k exp(i2ωt) + c.c.) and plugging this term into the SBEs, one obtained the complex contribution to the inversion as The second contribution δN B k enters due to intraband dynamics or Boltzman type transport.According to Eq. ( 18), the intraband dynamics can be incorporated by displacing the density matrix in k-space by ∆k.As a result of two subsequent translations conducted by two conjugate fields with the frequency of ω and −ω, the steady-state population difference experiences a minor change denoted by δN B k .By Taylor expanding the population, one obtains zero around the Dirac point, which confirms the slow relaxation dynamics suggested above.

Figure 1 .
Figure 1.(a) Low-energy band structure of graphene.The angle φk and the magnitude of the Bloch wavenumber k are defined in polar coordinate system in reciprocal space.(b)& (c) k-dependent population and coherence relaxation coefficients shown by γ (1) k and

Figure 2 .
Figure 2. Change in the conductivity of graphene seen by the probe field (with frequency ωp) in a pump-probe experiment for the pump frequencies (a) ωc = 80meV and (b) ωc = 800meV.The conductivity is normalized to σ 0 = e 2 /4 .Corresponding changes of the absorption coefficient of graphene for different intensities of the pump field are shown to the right of the conductivity plots.

Figure 3 .
Figure 3. (a) In the calculation of the semiperturbative nonlinear coefficient, the saturated region displayed by the dark blue disc is excluded (b-d) Kerr-type nonlinearity of graphene obtained from the analytical nonperturbative approach (yellow shaded curves) and the semiperturbative approach (i.e., |σ (3) (ω, ω, −ω)|, black dotted curves) plotted for different field magnitudes (i.e.E).The z-axis corresponds to |σ (3) (ω, ω, −ω)|/σ 0 where σ 0 = e 2 /4 .The two distinctive blue shaded regions are the saturation regions in E − ω plane due to, respectively, zero-detuning (light blue) and the strong interband coupling in the vicinity of the Dirac points (dark blue).Results for different Fermi levels (E f = 100, 200, 400meV) are displayed in (b-d).For low Fermi levels the nonlinear optical response is noticeably field-dependent.if the photon energy lies on the dark-blue disk semiperturbative coefficients would not exist and thus the black dotted curves are not extended within the low-energy saturation domain.In the zero-detuning saturation domain (light blue area), the semiperturbative analysis fails and it cannot follow the non-perturbative results.

Figure 4 .
Figure 4. Schematic of the effective-carrier dispersion in graphene.The vector k is the Bloch wavenumber with respect to the Dirac point.

Figure 5 .
Figure 5. Schematic representation of different three-photon processes contributing to the third-order nonlinear optics of graphene.The intraband type of dynamics is depicted by displacement of the distribution and the interband dynamics are shown by two-level transitions.

( 3 )
k .Assuming that the first photon ω p provides time variation for v

( 1 )
k (ω p ), the first nonzero oscillatory component of w k as well as the third harmonic of v(3)

Figure 6 .
Figure 6.Extraction of the microscopic relaxation coefficients for the optical excitation of ω = 80meV at (a) vF k = 10meV(low energy) and (b) for vF k = 100meV (high energy) points.

( 1 )Figure 7 .
Figure 7. Extraction of the microscopic relaxation coefficients for the optical excitation of ω = 800meV at (a) vF k = 50meV(low energy) and (b) for vF k = 500meV (high energy) points.

Figure 9 .
Figure 9. Schematic of the two-frequency pump-probe experiment.

Figure 10 .
Figure 10.change in the diagonal element of the conductivity at the probe frequency ωp (δσ(ωp)) due to the presence of the pump field at the frequency of (a) ωc = 80meV and (b) ωc = 800meV.The conductivity is normalized to σ0 = e 2 /4 .Variations of the absorption coefficient of graphene (normalized to its intrinsic value) versus the pump intensity for two frequencies of the probe field are shown in the insets of the figures.
Work at Waterloo has been supported by the Natural Science and Engineering Research Council of Canada (NSERC) and the Canada First Research Excellence Fund.Work at Chalmers has been supported by the Rune Bernhardsson Graphene fund, by Vetenskapsrådet under grant No. 2016-03603, and by the European Union's Horizon 2020 research and innovation programme under grant agreement No. 785219.