{'fraction_non_alphanumeric': 0.039375618908636056, 'fraction_numerical': 0.013159952050867774, 'mean_word_length': 4.328381005276312, 'pattern_counts': {'":': 0, '<': 2, '<?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': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comprehensive analysis of the optical Kerr coefficient of graphene

We present a comprehensive analysis of the the nonlinear optical Kerr effect in graphene. We directly solve the S-matrix element to calculate the absorption rate, utilizing the Volkov-Keldysh- type crystal wave functions. We then convert to the nonlinear refractive index coefficients through the Kramers-Kronig relation. In this formalism, the source of Kerr nonlinearity is the interplay of optical fields that cooperatively drive the transition from valence to conduction band. This formalism makes it possible to identify and compute the rates of distinct nonlinear processes that contribute to the Kerr nonlinear refractive index coefficient. The four identified mechanisms are two photon absorption, Raman transition, self coupling, and quadratic AC Stark effect. We also present a comparison of our theory with recent experimental and theoretical results.


I. INTRODUCTION
The Kerr optical effect describes the change in the refractive index in a material due to the presence of a strong optical field. It causes various nonlinear phenomena including self-focusing, self-phase modulation, and modulation instability. The Kerr effect can be used for optical bistability and all optical switching in the context of cavity nonlinear optics, which are the prominent building blocks for integrated and quantum photonics [1][2][3].
Graphene's Dirac cones exhibit a linear dispersion, making it an excellent material for variety of electrical applications (see the review articles [4,5]). More recently the "giant" optical nonlinearity of graphene has been explored by theory and experiments [6][7][8][9][10][11][12][13] (see a comprehensive review of recent progress in Cheng et al. [9]). The third-harmonic generation and the Kerr-type intensitydependent refractive index change have been investigated in those previous results. In general, near-band four-wave mixing and Kerr optical nonlinearity are similar processes while the third-harmonic generation is somewhat distinct process that involves waves that are widely separated in frequencies. Although there are a number of heuristic approaches including the Miller rule [14] to relate those two processes, the relation of the third-order susceptibilities describing the two processes are not rigorously proven, and this treatment generally assumes a direct-gap band structure with ω E g [15]. For graphene's optical nonlinearity, rigorous theories describing the third harmonic generation have been reported, including Mikhailov [16] who calculated the off-diagonal current density matrix for the interband transition and Cheng et al. [9], who used the semiconductor Bloch equation also to calculate the current density matrix. Although Cheng et al. briefly showed the Kerr nonlinear coefficient result based on twophoton absorption, it is known that other mechanisms also contribute to the Kerr nonlinear coefficient [17]. To understand the optical Kerr nonlinearity in graphene, it is thus necessary to exhaustively search for the physical mechanisms that contribute to the Kerr nonlinear coeffi-cient.
In this paper, we take the Sheik-Bahae approach that relies on the Kramers-Kronig relation to convert the absorption coefficients to dispersive coefficients, directly calculating the S-matrix element describing the transition from the valence to the conduction band, using the Volkov-Keldysh-type dressed-state wavefunctions [17]. The advantage of this method, compared to others such as calculating the current density operators, is that it allows for complete identification of the contributing physical mechanisms through a closed-form analytic expansion that includes an arbitrary level of wave mixing processes. Using this method, we identified a full set of physical processes that include various separate twophoton transitions as well as energy band corrections, which also naturally appear in the formalism. For the first time that we are aware of, we successfully calculated the Kerr nonlinear coefficient of graphene as a function of optical frequency and produced an analytical formula.
We first explain the strategy to calculate the Kerr nonlinear coefficient (section II). Once the formulation is laid out, we perform an exhaustive search for the two-photon transitions as well as energy band corrections (section III), producing a table of contributions. We then convert the contributions to finally calculate the Kerr nonlinear coefficients (section IV). A comparison between our theory and previous experimental and theoretical results is presented (section V), and conclusions follow (section VI).

II. CALCULATION METHOD
Our method to calculate the Kerr nonlinear coefficient utilizes the Kramers-Kronig relation where n is the refractive index, c is the speed of light, α r is the absorption coefficient in unit (m −1 ). The sym-arXiv:1605.03673v2 [cond-mat.mes-hall] 4 Aug 2016 bol P indicates that the integral is a Cauchy principal value. We note that the Kramers-Kronig theory is intrinsically a linear relation, allowing each contribution to be calculated separately and added up later. We follow Sheik-Bahae's argument [17] to transform the above into a perturbative calculation relating high order terms between the absorption coefficient and the refractive index. Suppose that the absorption coefficient changes due to a presence of perturbing source at frequency Ω. Then, the amount of change in absorption at frequency ω can be written as ∆α r (ω ; Ω). Eq. (1) implies that the induced change in the real refractive index ∆n(ω; Ω) can be still similarly linearly related to the change ∆α r (ω ; Ω): We can also calculate the self-refractive index change, setting Ω = ω: which describes the change of refractive index at frequency ω, caused by the light of the same frequency. If one adopts a perturbative approach, this equation implies that the perturbative expansion of the refractive index is related to that of the absorption coefficient, order by order. The Kerr nonlinear coefficient n 2 is defined as ∆n(ω) = n 2 (ω)I ω where I ω is the intensity of frequency ω. Hence, our task is to find an appropriate formula for ∆α r (ω , ω) ∝ I ω . Consider the amount of absorbed light at ω from the absorption by the graphene. The light intensity difference due to the absorption is δI ω = −αI ω , where α is the absorption rate. The absorption may change due to the presence of light at frequency ω. The change in absorption is ∆δI ω = −∆α(ω ; ω)I ω . Here, ∆α(ω ; ω) is the absorption rate change at ω , affected by the presence of light at ω. The generated number of freecarriers from absorption is N f = −A · (δI ω ) · dt/( ω ) where A is the area of graphene that is exposed to the light field, dt is the time of exposure. The change of the number of free-carriers due to the presence of light at ω is dN f = −A·(∆δI ω )·dt/( ω ). Then, the change in the rate of free-carrier density due to the presence of light at ω is, using dn f = dN f /A: We express such change perturbatively with respect to I ω , the intensity of light at ω, and take the leading term. Then, we can express the changed amount of the free-carrier density generation rate by (dn f /dt) = A(ω , ω)I ω I ω with an appropriate parameter A(ω , ω). Then, we can find the effective absorption coefficient, dividing ∆α(ω ; ω) by an effective atomic thickness d 0 of the graphene plane [9] to obtain Finally this will allow us to obtain the Kerr coefficient n 2 using the relation in Eq. (3): Hence, to calculate n 2 it is now apparent that we need to calculate the change in the free-carrier density generation rate, expressed as a bilinear term I ω I ω .
The free-carrier density generation rate is related to the transition rate of the electron state between the valence band state and the conduction band state via the interaction Hamiltonian H int = (e/m)p · A with electron charge e, mass m, the canonical momentum p, and the vector potential A of an external field: where W is the transition rate per area, D = 4 is the degeneracy factor of graphene (factor 2 from two disparate Bravais lattices and factor 2 from the two spin states), and W i→j (k, k ) is the transition rate from a Bloch wave with wavevector k in the band i to another with k in the band j. Here, f (i, k) represents the Fermi distribution (occupation probability) of band i, for a given temperature and doping. If we restrict ourselves to the case of an ultrafast optical pulse so that we can ignore the thermal recombination rate, we can approximately set W v→c (k, k ) W c→v (k , k), which leads to The transition rate W v→c (k, k ) is calculated through the S-matrix formalism [18]: where dt is the incremental time of exposure, which is identical with the previous dt. We use the first-order perturbation with respect to the interaction Hamiltonian and represent the S-matrix element: where t is the total time of exposure, |ψ c,v (k) are the Dirac ket of the conduction and the valence band states with the Bloch wavevector k, andH int (τ ) is the interaction Hamiltonian in the interaction picture from interaction Hamiltonian H int = (e/m)p · A(τ ).
A typical approach for calculating the third-order optical nonlinearity is to utilize the second order perturbation theory in the S-matrix formalism so that squaring it generates the third-order with respect to the interaction Hamiltonian. An alternative is to consider only the first-order perturbation and utilize the dressed state wavefunction at time τ [17]. For this, we need to obtain an expression for ψ c,v (k, r, τ ) = r|ψ c,v (k, τ ) .
Volkov first solved the Dirac equation for a free electron under the influence of an external field [19], relying on the fact that the lack of binding potential makes the interaction Hamiltonian commute with the free unperturbed Hamiltonian p 2 /2m. This allows for a relatively straightforward full solution. In a crystal lattice, however, one cannot use the same approach due to the presence of the binding potential of the lattice atoms. Keldysh in [20] proposed a hybrid Volkov-type wave function for a Bloch wave inside a crystal lattice. Concretely, Keldysh's modified Volkov-type wavefunctions are where u j (k, r) is the Bloch function that has the same period as the lattice, j = c, v is the band index for either conduction or valance band respectively, and Keldysh's wavefunction is obtained by time-integrating the Schrödinger equation after replacing the unperturbed Hamiltonian with energy eigenvalues (band energies). At the same time, one replaces the canonical momentum p with the crystal momentum k, which are not identical in general, thereby, making the interaction Hamiltonian commute with it. Obviously if t = 0, the wavefunction reduces to a usual Bloch wavefunction. If the external field A is zero, the time evolution of the wavefunction follows the phase rotation based on the band energies.
Here, replacing the mass with the effective masses m c,v is critical since the binding potential is effectively absorbed in the band-structure dependent effective mass as well as the band energies. As a result, the modified wavefunction appears as a free propagating wave with a modified energy through the first-order perturbation theory, periodically spatially modulated through the Bloch function u i (k, r). This is exactly what Bloch wavefunction is supposed to mean physically. This dressed state formalism, as Sheik-Bahae et al. pointed out, allows for a simple integral calculation of the S-matrix element, leading to identification of a complete set of physical mechanisms affecting the Kerr nonlinear coefficient, using only the first-order perturbation theory [17]. This is the beauty of Keldysh's treatment of dressed state of the lattice crystal Bloch wavefunctions. The result was extremely successful that the experimental results from many semiconductors including both wide band-gaps and narrow band-gaps, spanning several orders of magnitude in n 2 , match the calculated Kerr coefficient n 2 reasonably well (within a factor of 5).
Graphene is radically different from the bulk semiconductors in multiple aspects, and thus, we cannot adopt Sheik-Bahae's approach directly. Beside the fact that graphene is a 2D material, Sheik-Bahae et al. treated only the case for the typical quadratically dispersive materials with a non-zero bandgap energy, as shown in Eq. (12). This was the critical step, which allowed Sheik-Bahae et al. to be able to calculate the Kramers-Kronig integration. Unfortunately graphene has a linear dispersion, the effective mass diverges, and the bandgap energy is zero.
To overcome the difficulty in treating graphene while retaining the spirit of Keldysh's crystal lattice full solution, we study whether we can replace the Hamiltonian usually appearing in the time integrated exponent with C-numbers. Let us investigate this possibility by setting where E j (k) is the band energy and |ψ j (k) has a projection on to the position basis: the Bloch wavefunction We put this into Eq. (10) and insert the closure relation Noticing that m e i(k−k )·Rm = N δ k,k where R m is the mth lattice point and N is the total number of atoms, and the S-matrix element is thus nonzero only when k = k (lattice momentum conservation), we obtain where To explain the above equation, we note that there are four terms appearing in S cv after putting the closure relation. We also note that H int ∝ |A| and all the terms in S cv are at least the first order in H int . Since we are looking for the terms proportional to I ω I ω after squaring S cv , and the intensity I ω is proportional to |A| 2 , the terms we are interested in S cv , before squaring, are at most first order in H int . However, the terms involving ψ c (k)|e i dτ Hint |ψ v (k) (= B) does produce only the second-order or higher in H int , before squaring S cv , because of the following reason. When perturbatively expanded, the leading orders of B and B * are both the first order in H int . Therefore, when multiplied by H int and other terms, the terms containing B and B * result in the second order in H int , before squaring S cv . Again, since all terms are at least the first order in H int , the terms containing B and B * will result in the third order in H int or higher when squared. Therefore, the terms containing B and B * are discarded in the above equation, and we are left with only one term.
This explicit S-matrix calculation implies that we may utilize the following wavefunction, which is correct in the first order perturbation theory of the interaction Hamiltonian: The remaining task is to evaluate the matrix elements of the interaction Hamiltonian in the basis Let us first note that where v ij is the matrix element of the velocity operator Here, we approximated the interaction Hamiltonian only up to the first order perturbation of the field A, which is necessary to obtain the terms I ω I ω . This is the critical aspect of our calculation method that makes the whole calculation procedure tractable. The velocity operator, around the K-points (Dirac cones) of the graphene, is obtained previously in [9], where a sketch of derivation was presented. We present a concrete derivation in the appendix A since the elements of the velocity operator play an essential role in our calculation. The results around the Dirac cones are with v vc = v * cv . Here, v F is the Fermi velocity equivalent to 3at nn /2 with a = |a| where a is the vector connecting the atom in the Bravais lattice A to the nearest atom in the Bravais lattice B (a = 0.139 nm), and t nn ≈ 2.7 eV is the hopping energy between the nearest atom sites [21]; z is the direction of light propagation, perpendicular to the graphene 2D plane, and q = k − K. We note that the velocity matrix is Hermitian.
Next, let us assume explicitly that the external field is represented as where N ω is the number of light waves at different frequencies andâ is the unit vector such thatâ ·ẑ = 0.
Here, we assumed that all light fields at various wavelengths have the same polarization since we are primarily interested in the Kerr nonlinearity. The interaction Hamiltonian is, then, Then, the explicit wavefunction obtained by Eq. (18) is Finally we obtained an explicit form of the wave functions consisting in all C-numbers. This is one of the key steps in our procedure.
To evaluate the S-matrix element, we further use the Bessel function expansion: After simplifying the time integral into Dirac delta functions, we finally obtain an explicit formula for the Smatrix element: where This is the most general form of the S-matrix using the first-order perturbation theory based on the dressed state formalism. Now that we have established our calculation procedures and obtained the elements necessary to carry them out, we want to emphasize that what we are searching for is the term proportional to I ω2 in calculating the absorption coefficient for frequency ω 1 (see Eq. (5)). Since the absorption coefficient is derived from the free-carrier transition rate, which is proportional to the square of the S-matrix element in Eq. (27), what we are searching for in the S-matrix elements are the terms that are proportional to A 1 A 2 where A i is the field amplitude at frequency ω i , since the amplitude of the field and the intensity are related as I ωi = (cε 0 ω 2 i /2)|A i | 2 where ε 0 is the vacuum permittivity. These terms eventually contribute to the Kerr nonlinear coefficient through the Kramers-Kronig relation in Eq. (3). Therefore, we are primarily interested in searching for terms proportional to A 1 A 2 in the S-matrix.
In summary, we will first calculate the S-matrix element to obtain the transition rate from the valence to the conduction band via the interaction Hamiltonian, which is translated into the carrier density change rate and relates to the absorption coefficient. Then, we will evaluate the Cauchy principal value integration using the Kramers-Kronig relation, which finally leads to the Kerr nonlinear coefficient n 2 .

III. ABSORPTION COEFFICIENTS
In this section, we will calculate the absorption coefficients for both linear and nonlinear processes. We first verify that our formalism produces the well known result of the graphene linear absorption. Then, we proceed to the two-photon absorption and energy corrections, which are the essential stepping stones towards the Kerr nonlinear coefficient n 2 .

A. Linear absorption
The linear absorption is the case where only one light wave is involved in the transition. Hence, we set N ω = 1 and n 1 = 0 in the S-matrix element in Eq. (27) to obtain We note that the second delta function in Eq. (27) disappears since all frequencies are positive in our formalism.
Using the typical treatment of square of the delta function [22]: we obtain the carrier density rate from the linear absorption using Eq. (8) and Eq. (9) We will replace the summation of k into integral. For this, we need to evaluate the following using the previous result in Eq. (21): Converting the sum into integral ( For an undoped graphene, we have f vc (E) = tanh(E/2k B T ) where k B is the Boltzmann constant and T is the temperature. Let us use the relation between the field amplitude A 1 and the intensity of the field I ω1 , given by I ω1 = (cε 0 ω 2 1 /2)A 2 1 . Using the relation in Eq. (4) for the linear case, we finally obtain the linear absorption rate This exactly matches the well-known broadband graphene linear absorption 2.3 % in the optical frequencies ( ω 1 k B T ) [4,23]. It is noteworthy that the effect of doping, and hence, the non-zero Fermi energy E F , is dramatically pronounced for low frequencies (long wavelengths), which exhibit a non-analytic behavior at low temperature. The treatment in [24] explicitly approximated the Fermi distribution functions at low temperatures and demonstrated such non-conventional behavior at longer wavelengths ( ω 1 2E F ) in the low temperature limit. Such effect, however, is relatively small at optical frequencies.

B. Two-photon absorption
The case corresponding to the two photon absorption is N ω = 2 (that is, two light frequencies are involved). If we use the Bessel function's Taylor series expansion to the lowest order: we expect that, in the S-matrix element Eq. (27), n j corresponds to the number of times A j is multiplied. Recall that we are searching for the terms in the S-matrix element, which are proportional to A 1 A 2 . These terms are obtained if in Eq. (27), for i = 1, we set n 1 = 0 and n 2 = 1 while, for i = 2, we set n 1 = 1 and n 2 = 0. The delta functions enforce that only three cases are possible, satisfying: Figure 1 shows the transitions that the above three cases represent. The first transition describes the case where ω 1 , ω 2 < ω cv , which is the typical two-photon absorption. The second and the third transitions describe the Raman transitions where one frequency is larger than the band energy difference. We will first calculate the absorption coefficient for the case where both frequencies are smaller than the band transition energy. The S-matrix element is To proceed further to calculate the change of the freecarrier density generation rate, we must evaluate the following, using Eq. (21), where we considered the fact that v ii is orthogonal to v cv while they both reside in the graphene 2D plane. Using this and after a long but elementary calculation, we obtain (39) Converting the field amplitudes into the intensities and using Eq. (4), we obtain the absorption coefficient The degenerate case where ω 1 = ω 2 is of particular interest. It immediately follows from the above equation that the degenerate two-photon absorption rate (which is the absorption coefficient multiplied by d 0 ) is Note that the two photon absorption is proportional to the inverse fourth power of the optical frequency ω. This result agrees precisely with the reference [9].

C. Raman transition absorption
As we indicated in the previous section, another cooperative two-photon absorption is possible through the second, and the third case in Eq. (36). For the case where ω 2 > ω cv (the second equation in Eq. (36)), we obtain a term in the S-matrix that is proportional to A 1 A 2 by setting N ω = 2 and n 1 = 0, n 2 = 1. We ignore the case where either n 1 or n 2 becomes negative since those terms do not produce the terms proportional to A 1 A 2 . Also for the other case where ω 1 > ω cv , we have the S-matrix element proportional to A 1 A 2 by setting N ω = 2 and n 1 = 1, n 2 = 0. Then, for these two cases, we obtain where i, j = 1, 2 with i = j. We need to be careful in calculating the absorption coefficient since we will be using Eq. (4) where I ω1 / ω 1 appears on the right hand side. The result is that We note that ω 2 > ω cv > ω 1 corresponds to ∆α RT,1 r whereas ω 1 > ω cv > ω 2 corresponds to ∆α RT,2 r . Both absorption coefficients are thus always positive.

D. Self-coupling
Another term in S-matrix element that is proportional to A 1 A 2 is through expanding the zeroth order Bessel function into in the linear absorption calculation in Eq. (29). That is, when calculating the linear absorption at ω 1 , the presence of light of frequency ω 2 may affect the absorption through the second-order expansion term in the above equation. Physically the process can be understood as a reduction in the oscillation strength at ω 1 by renormalizing the interband coupling due to the acceleration of the free-carriers caused by the field at ω 2 [17]. The result is the modified S-matrix element from the linear absorption Eq. (29) Squaring it produces the term that is proportional to A 2 1 A 2 2 , while we discard the higher order term proportional to A 2 1 A 4 2 . Through the similar calculation procedure as other mechanisms, we obtain the absorption coefficient from the self-coupling (46) As expected, the self-coupling results in the reduction of the absorption (the quantity above is negative).

E. Quadratic AC Stark effect
The last remaining contribution that produces a term proportional to I ω1 I ω2 can be traced by recalling Eq. (33) where the argument inside the delta function plays a role. That is, integrating E δ(2E − 2E ) with respect to E = v F k results in a factor E . This implies that, if for some reason the argument of the delta function is modified to 2E − 2E + ∆E , the factor appearing after integration must be (E − ∆E /2), in the change of the free-carrier density generation rate as well as in the absorption coefficient. Particularly if such modification involves a term that is proportional to A 2 2 of light at frequency ω 2 , a new term that is proportional to A 2 1 A 2 2 appears, which must be counted for to calculate the Kerr nonlinear coefficient. The AC Stark effect is such a quadratic energy modification for both the conduction and the valence band. The energy shift adds more on the band energy difference by [25] where Ω is the Rabi oscillation frequency and ∆ is the detuning of ω 2 light from the main resonance at ω 1 .
Note that the band energy may increase or decrease, depending on whether one has blue or red detuning. This quadratic AC Stark shift produces a contribution that is proportional to A 2 1 A 2 2 in the change of the free-carrier density generation rate. The result appearing in the absorption coefficient is, after integral calculations, Note that the overall sign depends on the sign of ω 1 − ω 2 . Hence, the quadratic AC Stark effect increases (decreases) the absorption rate if ω 1 > ω 2 (ω 2 > ω 1 ), respectively.

F. Summary of absorption coefficients
Let us define a lumped parameter where we used d 0 = 0.33 nm [9]. We can express the various absorption coefficients using the following convention: (50) Table I summarizes our calculated absorption coefficients from the various two-photon mechanisms. It is worth noting that the contribution from the self-coupling and the quadratic AC Stark shifts are the corrections to the linear absorption calculated in Eq. (34), and thus the negative sign does not mean creation of photons. Also, the occupation probabilityf vc reflects the energy difference of the actual conduction and the valence bands involved in each process.

IV. KERR NONLINEAR COEFFICIENTS
We are now ready to calculate the Kerr nonlinear coefficients n 2 using the absorption coefficients in Eq. (50) with the TABLE I, through evaluating the Kramers-Kronig relation in Eq. (3). We need to carry out the Cauchy principal value integration with respect to ω 1 . Using the relation ∆n(ω 2 ) = n 2 I ω2 in Eq. (3) and using Eq. (50), we obtain, replacing ω 1 /t nn → x 1 and ω 2 /t nn → x 2 , We first calculate the contribution from the twophoton absorption. The presence of the occupation functionf vc complicates the calculation. In principle, one can evaluate using a full numerical integration. However, with a reasonable assumption, we can carry out the integration analytically. For undoped graphene, the occupation function is f vc (E) = tanh(E/2k B T ). Particularly when ω 1,2 are optical frequencies where ω 1,2 k B T with a reasonable temperature range well below the Fermi temperature, one can safely approximate f vc ( (ω 1 + ω 2 )/2) 1. Then we obtain .
(53) One can express the Kramers-Kronig relation between the real and the imaginary values of the susceptibility χ = χ + iχ : where the real part χ is related to ∆n(ω) and the imaginary part χ is related to ∆α r (ν; ω). This is the same expression as Eq. (1) and Eq. (3). Due to the condition that χ (−ω) = −χ (ω) that comes from the realness of the susceptibility χ(t) (i.e., χ * (ω) = χ(−ω)), any singularity occuring around ν = 0 is canceled due to the antisymmetry between ν = 0 + and ν = 0 − of the integrand. Hence, on these physical grounds we can safely discard any mathematical singularity at x 1 = 0 in the evaluation of the Cauchy principal integral Eq. (53). Integrating with respect to x 1 in the separate intervals [x 0 , x 2 − ] and [x 2 + , ∞), taking the limit x 0 → 0 and then → 0, and discarding the singularity occurring at x 1 = 0, we calculate From this we obtain the contribution of the two-photon absorption to the Kerr nonlinear coefficient as We note that this quantity is positive if ω 2 < t nn , but negative otherwise. We also note that the unit of this quantity is (m 2 /W) such that n 2 (ω 2 )I ω2 becomes unitless.
Calculating the Kerr coefficient from the two Raman transition processes requires to evaluate the following: Since we are concerned about optical frequencies ω 1 , ω 2 , it is also a reasonable assumptions to set where θ(x) is the Heaviside step function. Then, the above reduces to .
Ignoring the singularity occurring at x 1 = 0, we obtain Hence, we obtain We note that this quantity is positive if ω 2 > 0.41t nn , and otherwise negative. To calculate the contribution from self-coupling, we need to evaluate (62) After approximating f vc (x) θ(x), it is easy to verify that the above integral vanishes.
Next, we evaluate the contribution from the quadratic AC Stark shift, which requires the evaluation of the following integral.
We first approximate f vc (x) θ(x). Next, we perform the Cauchy principal value integral with respect to x 1 in the two separate intervals [x 0 , x 2 − ] and [x 2 + , ∞). Adding the results together, discarding the singularity occuring at x 1 = 0, and also sending x 0 → 0 and then → 0 produces a converging analytic result, except for one term involving 1/(x 2 2 ). The origin of this term is the divergence of the quadratic AC Stark shift near x 2 ≈ x 1 in Eq. (48), which is unphysical since the form we used is valid only when the detuning is large. With a more physically accurate formula involving the natural line broadening of the upper state, it is well known that the AC Stark shift for the case x 2 = x 1 is zero. Moreover, the contribution from x 1 = x − 2 is exactly canceled out by the contribution from x 1 = x + 2 . Hence, in the Cauchy principal value integral we discard this diverging term. The obtained result is Let us define a lumped parameter We can represent the various Kerr nonlinear coefficients using the following convention: (ω 2 ) = β 2 t nn ω 2 5 2 ln 2 + 6 ln t nn ω 2 .
We note that n total 2 changes sign at ω 2 = 1.26t nn (equivalently at 365 nm). If ω 2 > 1.26t nn , the Kerr coefficient is positive and, otherwise, negative.  comparison of the magnitudes of the contributing mechanisms to the Kerr nonlinear coefficient. At wavelengths longer than 0.6 µm, the contribution from the two photon absorption is the largest. However, contributions from the Raman transitions and the quadratic AC stark effect are not small. It is clearly shown that the Kerr nonlinear coefficient is not singly determined by the two-photon absorption. A comprehensive calculation is thus necessary.

V. DISCUSSION
We first discuss the valid region of our calculations. Recall that we have assumed the linearized dispersion around the Dirac cones, which is explicitly used in de- riving the analytic velocity matrix element in Eq. (21). The Kramers-Kronig relation in Eq. (3) integrates ω 1 ∈ [0, ∞). FIG 5 in the appendix clearly shows that the band transition energy ω 1 larger than 4 eV (equivalently, wavelengths shorter than 310 nm) deviates from the linear dispersion. Hence, one must resort to numerical integration of the Cauchy principal integral, in order to address the wavelengths below 310 nm. Additionally, one clear consequence of our linearized assumption used in the Cauchy principal value integral is that the Kerr nonlinear coefficient from the self-coupling became zero. If we take the same approach to calculate the linear refractive index n 0 through the Kramers-Kronig relation Eq. (1), utilizing the linear absorption in Eq. (34), it is straightforward to verify that we obtain n 0 = 1. However, a known experimental result exhibited that the refractive index of the graphene approaches that of the graphite n 0 ≈ 2.6 [5]. Hence, it is expected that our result may be different from the true values by a factor of unity magnitude. To compensate for this error, one may attempt to introduce a distortion function inside the Cauchy principal value integration for the high frequency components and artificially impose the consistency with the experimental result. However, this approach is somewhat discouraged since it is difficult to make an explicit formula for such a distortion function and carry out consistently all the other Cauchy principal value integrations for higher order mechanisms. Furthermore, our formalism did not include the participation of other higher or lower energy bands that may have small effects. Therefore, it should be clearly understood that our analytical final result in Eq. (67) provides rather an order of magnitude estimate.
We compare our theoretical result with previously reported results. We carefully selected the results that are only relevant to the Kerr nonlinear coefficient n 2 , particularly excluding the results based on the third- (m 2 /W)) using four-wave mixing technique [7]. However, Cheng et al. [9] pointed out that the calibration of this result was not correct and recalculated based on the correct calibration, bringing down Hendry's result to 1 × 10 −17 (m 2 /W), which agrees with our result in an order of magnitude estimate. Zhang et al. reported n 2 10 −11 (m 2 /W) using Z-scan method at 1.55 µm [8], which is three orders of magnitude larger than our theory. This experiment was performed on a graphene sitting on a quartz substrate, and it was hard to distinguish the pure graphene response from the combined response of the graphene, the interface, and the substrate compound. On the other hand, the result of Gu et al. in [26] was obtained from a hybrid waveguide consisting of silicon and a mono-layer graphene. Utilizing the four-wave mixing spectra, Gu et al. measured the combined Kerr nonlinear coefficient of the entire waveguide as n 2 = 4.8 × 10 −17 (m 2 /W) at 1.56 µm. Utilizing their calculation method based on the field filling-factors, one can reverse-estimate the sole contribution of the monolayer graphene. The resulting estimate of a mono-layer graphene is n 2 = 6.2 × 10 −14 (m 2 /W), which agrees with our theoretical result within an order of magnitude. For theoretical results, Cheng et al. calculated n 2 ∼ 10 −15 (m 2 /W) at 1.56 µm as an order of magnitude estimate, which agrees well with our result. We note that Cheng et al. considered only the two-photon absorption mechanism omitting the Raman transition and the quadratic AC Stark shift.
(70) FIG 4 shows the plot of R(ω). At shorter wavelengths, the incoherent two-photon absorption clearly dominates whereas at longer wavelengths the coherent transition picks up and become significant.

VI. CONCLUSION
We reported a comprehensive theoretical analysis on the Kerr nonlinear coefficients of graphene at optical wavelengths. Our method is based on the Volkov-Keldysh-type dressed state treatment of the lattice wavefunctions, which allowed the first order perturbation to describe all relevant mechanisms that contribute to the Kerr nonlinear coefficient, namely, the two-photon absorption, the Raman two-photon transition, the selfcoupling, and the quadratic AC Stark shift effect. We obtained a neat analytical formula for n 2 . If desired, one can obtain more accurate result by numerical calculation using the full-band model, based on our formalism. Nevertheless, our analytical solution provides an order-of-magnitude estimate of the Kerr nonlinear coefficient. When compared with previous results, our theory predicts around the average order-of-magnitudes of the reported results, partly due to the widely disagreeing experimental results caused by the complexity of the measurement setup. We also analyzed the ratio between the coherent and the incoherent processes using the optical nonlinearity of the graphene, which suggests that long wavelengths in deep infrared or THz waves are much more preferred to perform cQED experiments using graphene.
where s = ±1 if j = c, v respectively, γ k = 1 + e ik·a1 + e ik·a2 is the structural factor with a 1 = √ 3a(( √ 3/2)x − (1/2)ŷ) and a 2 = √ 3a(( √ 3/2)x + (1/2)ŷ) are the primitive lattice vectors, a = 0.139 nm is the distance between the nearest neighbor atomic sites, and φ i (r) is the |2p z orbit wavefunction of a single atom centered at the origin. The energy eigenvalues are E j = s|γ k |t nn with the hopping energy between the nearest neighbor sites t nn = 2.7 eV [21]. Figure 5 shows the energy band structure of the graphene around one of the K-points. Near the K-point, the dispersion is circularly symmetric. Moreover, approximately up until the energy difference E c − E i becomes 4 eV, the dispersion is quite linear.
Threfore, we can linearize, around the K-points, the structural factor to be γ k ≈ ( v F /t nn )(iq x + q y ) where (q x , q y ) = q = k − K. Here v F is the Fermi velocity given by v F = 3at nn 2 . (A2) The energy eigenvalues of the conduction and the valence bands become We first resolve the off-diagonal velocity matrix element v cv = ψ c |v|ψ v . Let us define the Berry phase ξ cv ≡ R0 d 2 r u c (k, r)∇ k u v (k, r).
Using the Bloch function in Eq. (A1), one can show that Using the basis {|Φ A , |Φ B } where r|Φ α = (1/ √ N ) m e ik·(Rm−r) φ j (r − R m − τ α ) with τ A = 0 and τ B = a is the vector from an atom in Bravais lattice A to the nearest atom in Bravais lattice B, we can express the kets where s = ±1 respectively for j = c, v. Using this representation, it is straightforward to calculate whereẑ is the direction of light propagation, which is perpendicular to the graphene plane. The last step is to calculate using Eq. (A3), Eq. (A5), and Eq. (A7): In addition, it is straightforward to verify that v vc = v * cv . For the diagonal terms, we utilize the following group velocity of the Bloch wave: where i = c, v. This immediately produces the result in Eq. (21).

ACKNOWLEDGMENTS
R. Hamerly is supported by a seed grant from the Precourt Institute for Energy at Stanford.
Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energys National Nuclear Security Administration under contract DE-AC04-94AL85000.