Nonlinear optical conductivity of a two-band crystal I
D. J. Passos, G. B. Ventura, J. M. B. Lopes dos Santos, J. M. Viana Parente Lopes
NNonlinear optical conductivity of a two-band crystal I.
D. J. Passos, G. B. Ventura, J. M. B. Lopes dos Santos, and J. M. Viana Parente Lopes
Centro de F´ısica das Universidades do Minho e Porto, Departamento de F´ısica e Astronomia, Faculdade deCiˆencias, Universidade do Porto, 4169-007 Porto, Portugal
The structure of the electronic nonlinear optical conductivity is elucidated in a detailed study ofthe time-reversal symmetric two-band model. The nonlinear conductivity is decomposed as a sum ofcontributions related with different regions of the First Brillouin Zone, defined by single or multiphotonresonances. All contributions are written in terms of the same integrals, which contain all informationspecific to the particular model under study. In this way, ready-to-use formulas are provided that reducethe often tedious calculations of the second and third order optical conductivity to the evaluation of asmall set of similar integrals. In the scenario where charge carriers are present prior to optical excitation,Fermi surface contributions must also be considered and are shown to have an universal frequencydependence, tunable by doping. General characteristics are made evident in this type of resonance-basedanalysis: the existence of step functions that determine the chemical potential dependence of electron-hole symmetric insulators; the determination of the imaginary part by Hilbert transforms, simpler thanthose of the nonlinear Kramers-Kr¨onig relations; the absence of Drude peaks in the diagonal elementsof the second order conductivity, among others. As examples, analytical expressions are derived for thenonlinear conductivities of some simple systems: a very basic model of direct gap semiconductors andthe Dirac fermions of monolayer graphene.
Early in the history of nonlinear optics, shortly after the discovery of second harmonic generation [1] and two-photon absorption [2], a general framework was set up to describe the variety of observed optical phenomena.This framework is based on the observation that even for the intense laser light required in these experiments,the electric fields involved are often much weaker than the atomic electric fields that bind the electrons,enabling a perturbative treatment of the light-matter interaction [1]. By means of an expansion of the electriccurrent (polarization) in powers of the optical fields, it is possible to define nonlinear optical conductivities(susceptibilities) and describe much of nonlinear optics [3]. These nonlinear response functions define whicheffects are present in a given medium, their frequency dependence and their efficiencies. For the purposes ofthis work, it will be assumed that the electronic contribution to the conductivity is dominant.Calculations of the nonlinear optical conductivity usually follow a semiclassical approach based on thedensity matrix formalism. The simplest treatments neglect electron-electron interactions and are withinthe electric dipole approximation; these assumptions are adopted here as well. For atomic and molecularsystems, general expressions for the conductivities are easily derived from perturbation theory and can, bythe use of inherent permutation symmetries , be reduced to a single term [3]. Despite the apparent simplicityof the nonlinear optical response in these systems, the requirement of knowledge of the exact eigenstates andenergies of the unperturbed Hamiltonian H , together with other complications such as local field effects,makes calculations for any but the simplest molecules less straightforward. Still, semi-empirical “sum-over-states” calculations have been successful in achieving quantitative agreement with experiment [5, 6].The treatment of the nonlinear optical response of a solid is, however, more complex. Integration overthe First Brillouin Zone (FBZ) is necessary and defining the perturbation matrix elements in the Blochbasis can be a nontrivial task. The first full band structure calculations for semiconductors were performedmore than a decade after the general expressions for the nonlinear conductivity were derived [4, 7–10]. A We are referring here to overall permutation symmetry as applied to σ βα ...α n (¯ ω , ..., ¯ ω n ) with complex frequencies. See [4]and pages 70-74 of [3] a r X i v : . [ c ond - m a t . o t h e r] F e b rst look at the minimal coupling Hamiltonian ˆ H (ˆ r , ˆ p + e A ( t )) used in early attempts to study nonlinearoptics in solids suggests that a crystal could be seen as a collection of independent atomic systems, onefor each k in the FBZ. This is indeed true when the entirety of the band structure is taken into account.However, working with an infinite number of bands is hardly practical and when we inevitably settle with anapproximate model with a finite number of bands, this description is shown to be inadequate and filled withunphysical infrared divergences [11, 12]. Only recently has this problem been resolved with a reformulationof the perturbation theory that extends its validity to systems with finite number of bands [13–16].In the nineties, the difficulties associated with the use of the minimal coupling Hamiltonian (also called“velocity gauge”) led to the development of an alternative, if equivalent, perturbation theory, correspondingto a different gauge choice (termed “length gauge”) in representing the optical field [17]. In this gauge, theHamiltonian has the form ˆ H = ˆ H + e ˆ r · E ( t ), where e is the electron charge, ˆ r is the position operatorand E is the classical optical field. Expressions were derived for the second order conductivity of coldsemiconductors that are applicable to models with a finite number of bands and no spurious divergences arefound in this formulation [18, 19]. The length gauge approach has since become a well established methodof obtaining nonlinear optical response functions [20–23].In the length gauge, both intraband and interband transitions must be taken into account and aretransparently expressed in the structure of the perturbation theory, which contains, as particular cases, thedynamics of atomic systems and of the free carriers single-band motion, but is more general, and complex,than either [17]. The complexity stems from expressing the perturbation in terms of a position operatorwhich takes the form of a covariant derivative in the Bloch representation [12, 17, 24]. The successiveapplication of derivatives as we move to higher orders in perturbation theory leads to unwieldy expressionsfor the nonlinear conductivity and lengthy calculations even for the simplest models, the only ones for whichanalytical calculations are even attempted. Despite this, the results sometimes show surprising simplicityand structure. As an example, the third order conductivity of the system of massless Dirac fermions foundin monolayer graphene has the form [22], σ xxxx ( ω, ω, ω ) = C ω (cid:18) − G (cid:18) (cid:126) ω | µ | (cid:19) + 64 G (cid:18) (cid:126) ω | µ | (cid:19) − G (cid:18) (cid:126) ω | µ | (cid:19)(cid:19) (1)with G ( x ) ≡ Θ( | x | −
1) + iπ log (cid:12)(cid:12)(cid:12)(cid:12) − x x (cid:12)(cid:12)(cid:12)(cid:12) (2)where ω is the frequency of the incident monochromatic optical field, µ is the chemical potential, v F is theFermi velocity and C is a constant: C = v F e / (cid:126) .In physics, when elaborate and extensive calculations are required to derive simple and elegant results,one is led to suspect that a simpler and more insightful way to express the theory exists. This is theperspective we take here. Eq. 1 has the third order conductivity broken up into pieces that are relevant indifferent regions of the spectrum, depending on whether one-, two- or three-photon frequencies are closerto matching the “effective gap” given by 2 | µ | . This suggests that a resonance-based decomposition of theconductivity might be possible in general, leading to a more direct derivation of analytical results and easierinterpretation of the underlying physics.In exploring this possibility, we confine ourselves to the study of the two-band model, the solid-stateanalogue of the two-level atom. In similar spirit to the theoretical investigations of the nonlinear optics oftwo-level atoms during the seventies [3, 25], we expect a study of the two-band model to provide a firmfoundation for later investigations of more general systems, to allow the central concepts to emerge moresimply and to have a wide range of applicability, encompassing any situation where the incident photonfrequencies connect a single pair of conduction ( c ) and valence ( v ) bands.As described in previous work [12, 13, 17], for the purpose of understanding nonlinear optical properties,the specification of the band structure (cid:15) k a and the non-abelian Berry connection A k ab , with a and b as bandindices running over { c, v } , defines the electronic system under study. For simplicity, it is assumed that thesystem is time-reversal invariant: (cid:15) − k a = (cid:15) k a and A − k ba = A k ab (for the right choice of phases in the Bloch In the relaxation-free limit. A k ab (an overly simplistic description of semiconductors) inSection 5 and for the system of massless Dirac fermions present in monolayer graphene in Section 6. Themain conclusions are summarized in Section 7. In a system of non-interacting electrons moving through the periodic potential of a crystal, the chargemotion that is perturbed by a passing light wave can be appropriately described, if the optical fields E ( t )are sufficiently weak, by performing an expansion of the current density J ( t ) in a power series, J ( t ) = J (1) ( t ) + J (2) ( t ) + · · · + J ( n ) ( t ) + . . . (3)The linear term in this expansion provides the usual definition of the optical conductivity, J β (1) ( t ) = (cid:90) + ∞−∞ σ βα ( t − t (cid:48) ) E α ( t (cid:48) ) dt (cid:48) (4)with an implicit summation over repeated tensor indices α . Similarly, we can define nonlinear conductivitiesby taking into account higher powers of E , J β ( n ) ( t ) = (cid:90) + ∞−∞ · · · (cid:90) + ∞−∞ σ βα ...α n ( t − t , . . . , t − t n ) E α ( t ) . . . E α n ( t n ) dt . . . dt n (5)This defines the constitutive relation J ( E ) in the time-domain. Since we shall be concerned with thepole structure of the nonlinear conductivity and resonances are more naturally discussed in the frequencydomain, we ought to rewrite the previous equation, J β ( n ) ( t ) = (cid:90) + ∞−∞ · · · (cid:90) + ∞−∞ dω π . . . dω n π σ βα ...α n ( ω , . . . , ω n ) E α ( ω ) . . . E α n ( ω n ) e − i ( ω + ··· + ω n ) t (6)with σ βα ...α n ( ω , . . . , ω n ) ≡ (cid:90) + ∞−∞ · · · (cid:90) + ∞−∞ σ βα ...α n ( t , . . . , t n ) e iω t . . . e iω n t n dt . . . dt n (7)Due to causality, the time-domain conductivity is nonzero only for positive times: σ ( . . . , t i , . . . ) = 0 if t i <
0, setting a lower limit in the integration range of Eq. 7 and an upper limit in Eqs. 4 and 5. This isseen explicitly in the expressions derived from perturbation theory (Eqs. 10 and 11).3 .1 Density matrix formalism
The dynamics of the current density, when thermally averaged, is given by the density matrix ˆ ρ , J ( t ) = Tr (cid:16) ˆ J ˆ ρ ( t ) (cid:17) (8)whose time evolution follows the Liouville equation, i (cid:126) ∂ t ˆ ρ = [ ˆ H, ˆ ρ ] = [ ˆ H , ˆ ρ ] + e [ˆ r α , ˆ ρ ] E α ( t ) (9)To perform the expansion in Eq. 3, the density matrix itself must be expanded as a powers series in theoptical fields. Perturbative solutions to Eq. 9 can then be found and replaced in Eq. 8. When solving Eq. 9,it is assumed that in the infinite past, when the perturbation was absent, the system was in an equilibriumdescribed by the Fermi-Dirac distribution ˆ ρ ( t = −∞ ) = ˆ ρ ≡ f ( ˆ H ).At order n , J β ( n ) ( t ) = (cid:18) − i e (cid:126) (cid:19) n (cid:90) t −∞ dt n · · · (cid:90) t −∞ dt Tr (cid:16) ˆ J β [ˆ r α n I ( t n − t ) , . . . [ˆ r α I ( t − t ) , ˆ ρ ] . . . ] (cid:17) E α ( t ) . . . E α n ( t n )(10)where the subscript I stands for the interaction representation r I ( t ) = e iH t/ (cid:126) r e − iH t/ (cid:126) .From Eqs. 5 and 10, we arrive at the nonlinear conductivity, σ βα ...α n ( t , . . . , t n ) = (cid:18) − i e (cid:126) (cid:19) n Tr (cid:16) ˆ J β [ˆ r α n I ( − t n ) , . . . [ˆ r α I ( − t ) , ˆ ρ ] . . . ] (cid:17) Θ( t − t ) . . . Θ( t n − − t n )Θ( t n )(11)As already mentioned, we are interested in the frequency domain version of Eq. 11. It could be derivedfrom Eq. 11 by a Fourier transform (Eq. 7) or by solving Eq. 9 directly in the frequency domain. The tracein Eq. 11 is to be evaluated over the eigenstates of the unperturbed Hamiltonian H , the Bloch functions ψ k a with energies (cid:15) k a ; the trace involves not only a sum over discrete band indices a and b but also an integrationover a continuum variable k that takes any value in the FBZ. The end result is that the nonlinear opticalconductivity takes the form, σ βα ...α n ( ω , . . . , ω n ) = e n (cid:90) d d k (2 π ) d (cid:88) a,b J β k ba (cid:126) ω + · · · + (cid:126) ω n − ∆ (cid:15) k ab × (cid:20) r α n , . . . (cid:126) ω + (cid:126) ω − ∆ (cid:15) k ◦ (cid:20) r α , (cid:126) ω − ∆ (cid:15) k ◦ [ r α , ρ ] (cid:21) . . . (cid:21) k ab (12)with ◦ as the Hadamard, or element-wise, product, d as the dimensionality of the system and ∆ (cid:15) k ab ≡ (cid:15) k a − (cid:15) k b . A more detailed derivation of Eq. 12 can be found elsewhere [12].Aside from the integration over the FBZ, the expression for the nonlinear conductivity in Eq. 12 is identicalto that of atomic systems. The essential differences emerge when attempting to write the position operatorin the basis of eigenstates of H . While in the case of atomic systems ˆ r has a well-behaved representation,in crystals the situation is analogous to that of a free particle, where a representation based on momentumeigenstates is used. In the momentum basis, the position operator takes the form of a derivative. In theBloch basis, ˆ r is the covariant derivative [12, 24], r α k ab = i D α k ab ≡ i (cid:18) δ ab ∂∂k α − i A α k ab (cid:19) (13)In their pioneering work [17], Aversa and Sipe noted that the representation of the position operatorshould not cause any real difficulties as the operator occurs only in well-defined commutators,[ˆ r α , ˆ O ] k ab = i [ ˆ D α , ˆ O ] k ab = i ∂ α O k ab + [ ˆ A α , ˆ O ] k ab (14)4here ˆ O stands for a generic observable with well-defined matrix elements, diagonal in k .In particular, the matrix elements of the current density, J β k ba = − e ( − i ) (cid:126) [ˆ r β , ˆ H ] k ba = − e (cid:126) [ ˆ D β , ˆ H ] k ba = − e (cid:126) (cid:16) δ ab ∂ β (cid:15) k a − i A β k ba ∆ (cid:15) k ab (cid:17) (15)with the use of the commutation relation (cid:2) r β , r α (cid:3) = 0 (Appendix A). A careful inspection of the general formula derived in the previous section raises some subtle issues. Forany frequency ω for which a resonance exists, the denominators in Eq. 12 are strictly undefined. Thisdifficulty traces back to Eq. 7 where the existence of a frequency domain nonlinear conductivity relies onthe convergence of the Fourier transform. When ignoring any kind of relaxation mechanism, the response toan impulse given at an instant of time can last forever: σ ( t = + ∞ ) (cid:54) = 0 and the Fourier transform diverges.This problem can be circumvented by extending the definition in Eq. 7 to complex frequencies in the upperhalf-plane [3], σ βα ...α n (¯ ω , . . . , ¯ ω n ) = (cid:90) + ∞−∞ · · · (cid:90) + ∞−∞ σ βα ...α n ( t , . . . , t n ) e i ¯ ω t . . . e i ¯ ω n t n dt . . . dt n = (cid:90) + ∞−∞ · · · (cid:90) + ∞−∞ (cid:0) σ βα ...α n ( t , . . . , t n ) e − γt . . . e − γt n (cid:1) e iω t . . . e iω n t n dt . . . dt n (16)with ω = Re(¯ ω ) and γ = Im(¯ ω ).The extension to complex frequencies can be interpreted in two ways. One is to consider the responsefunction in Eq. 16 to be associated not to monochromatic waves, but to fields that are adiabatically switchedon in the infinite past E ( ω ) e − iωt e γt . A different perspective is to look at complex frequencies as a simplephenomenological method to introduce relaxation into the system. As stated in Eq. 16, the nonlinear con-ductivity in the (real) frequency domain σ (cid:48) ( ω , . . . , ω n ) ≡ σ (¯ ω , . . . , ¯ ω n ) is obtained from a Fourier transformof a time domain response function that has the form σ (cid:48) ( t , . . . , t n ) = σ ( t , . . . , t n ) e − γt . . . e − γt n and satis-fies the condition σ (cid:48) ( t i = ∞ ) = 0. This approach to relaxation is most certainly too simplistic to properlyaccount for all the possible relaxation mechanisms that are observed in experiments, but it provides a directand easy way to obtain sensible answers and it has advantages relative to the traditional approach of addingan phenomenological term to the equation of motion (Eq. 9) [13, 26]. For simplicity, we here take the imag-inary part of the frequencies to be a constant γ but more generally we could have γ = γ ( ω ). It would onlybe required that the function be even, in order for the reality condition to be maintained.When the relaxation-free limit is considered, where the imaginary parts of the frequencies are taken tozero from above, the integrand in Eq. 12 can always be defined as a distribution by making use of theSokhotski-Plemelj theorem, (cid:90) d d k (2 π ) d g k (cid:126) ¯ ω − ∆ (cid:15) k γ → + −−−−→ − (cid:90) d d k (2 π ) d g k (cid:126) ω − ∆ (cid:15) k − i π (cid:90) d d k (2 π ) d g k δ ( (cid:126) ω − ∆ (cid:15) k ) (17)For an atomic system, the distribution will be defined relative to an integral over the frequencies (Eq. 6)and relies on the condition that the spectral width of E ( ω ) be much greater than γ . For a crystal, thedistribution is accommodated by the presence of an integration over the FBZ and no restrictions must beapplied to the optical fields considered. By taking the limit γ → + in the expression σ βα ...α n (¯ ω , . . . , ¯ ω n ) = e n (cid:90) d d k (2 π ) d (cid:88) a,b J β k ba (cid:126) ¯ ω + · · · + (cid:126) ¯ ω n − ∆ (cid:15) k ab (cid:20) r α n , . . . (cid:126) ¯ ω − ∆ (cid:15) k ◦ [ r α , ρ ] . . . (cid:21) k ab (18)the nonlinear conductivity can be defined as a regular function for real frequencies and vanishing relaxation.It is implicitly assumed in this reasoning that no more than a single denominator diverges at a given k (seeSection 4). 5n Eq. 18, the nonlinear conductivity can be further extended into the lower half-plane by analyticcontinuation [4]. In this way, Eq. 18 provides a valid expression over the entire complex frequency plane,even in regions where Eq. 7 no longer applies and the response function is, therefore, no longer physical. From Eq. 18, a more explicit form of the nonlinear conductivity can be derived by expanding out all thecommutators and performing all the required differentiations (that follow from Eq. 14), resulting in a lengthyand rather cumbersome formula. This is the usual starting point when computing the nonlinear opticalresponse functions of semiconductors and other materials. Numerical integration is necessary except forsome cases where a low-energy description exists with very simple dispersion relation and eigenstates. Forthese systems, analytical calculations are sometimes possible but still often rather long and complicated. Inthis section, we attempt to bring some simplicity and clarity to the structure of the nonlinear conductivity,by separating out terms whose resonances are located in different regions of the FBZ.Since we are restricting ourselves to the analysis of a two-band system, there is a single (nonzero) energydifference in the denominators of Eq. 18, ∆ (cid:15) ab = ± ∆ (cid:15) cv , allowing for a partial fraction decomposition intoterms with a single denominator to be integrated, ( (cid:126) ¯ ω + · · · + (cid:126) ¯ ω i ± ∆ (cid:15) cv ) − with i ∈ { , . . . , n } . These termswe denote by σ βα ...α n i (¯ ω , . . . , ¯ ω i , . . . , ¯ ω n ) as they are associated with resonances involving an i number ofphotons . We shall see later how the real part of these contributions is entirely described by the propertiesof the crystal in the vicinity of regions of the FBZ where the resonance condition (cid:126) ω + · · · + (cid:126) ω i − ∆ (cid:15) cv = 0is met. Some terms will involve poles of higher orders, but these can reduced back to simple poles with anintegration by parts or, equivalently, by making use of the identities in Appendix B. Finally, there will beterms where the application of the position operator resulted in derivatives of the Fermi-Dirac distribution.These terms will be treated separately and are denoted by σ βα ...α n F (¯ ω , . . . , ¯ ω n ). An explicit application ofthe procedure outlined here can be found in Appendix C, where the second order conductivity is treated indetail. For now, we state it generally, σ βα ...α n (¯ ω , . . . , ¯ ω n ) = σ βα ...α n F (¯ ω , . . . , ¯ ω n ) + σ βα ...α n (¯ ω , . . . , ¯ ω n ) + · · · + σ βα ...α n n (¯ ω , . . . , ¯ ω n ) (19)The various pieces of Eq. 19 will be made explicit in the following sections, but it is useful to first inspecttheir structure in general terms. The one-photon contribution can always be written as σ α (¯ ω , · · · , ¯ ω n ) = (cid:88) j (cid:88) p C p j (¯ ω , · · · , ¯ ω n ) Π p ( α ) j (¯ ω ) (20)where all tensor indices where condensed into one α ≡ βα . . . α n and p stands for permutation. The sum in p implies p ( α ) runs over all permutations of α , with a specific coefficient for each permutation applied. Thecoefficients C p j (¯ ω , · · · , ¯ ω n ) are specified in the following sections for the linear, second order and third orderconductivities ( n = 1, 2 and 3, respectively), where it is observed that most of these coefficients are zero,making only a small number of permutations necessary in practice. The coefficients are independent of thedetails of the system under consideration (they depend solely on the optical frequencies). All dependenceon material properties in the sum of Eq. 20 is in the integrals Π αj that take the general form,Π αj (¯ ω ) = (cid:90) d d k (2 π ) d (cid:88) a,b g αj ( A , (cid:15) ) ab (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (21)with ∆ f ba ≡ f ( (cid:15) b ) − f ( (cid:15) a ) as a difference in Fermi functions and g αj as a set of functions, labeled by j = 1 , , . . . , that depend on the energies and their derivatives and on the non-abelian Berry connection A and its generalized derivatives. For notational ease, the k label has been dropped.Similarly, for the two-photon contributions, In due rigor, we work with a classical electromagnetic field and there are no photons present. It should however be clear thatwhen a proper quantum treatment is made (cid:126) ω i is the energy of an incident photon, thereby justifying the used nomenclature. α (¯ ω , · · · , ¯ ω n ) = (cid:88) j (cid:88) p C p j (¯ ω , · · · , ¯ ω n ) Π p ( α ) j (¯ ω + ¯ ω ) (22)and the generalization is obvious at this point, σ αi (¯ ω , · · · , ¯ ω n ) = (cid:88) j (cid:88) p C pij (¯ ω , · · · , ¯ ω n ) Π p ( α ) j (¯ ω + · · · + ¯ ω i ) (23)Since all contributions involve a combination of the same integrals with a changing argument, the calcu-lation of the nonlinear conductivity is reduced to the evaluation of the integrals in Eq. 21. The complexityand number of integrals to be evaluated increases with the order n of the nonlinear conductivity, but theyalways retain the general form of Eq. 21 with varying g functions. For finite γ , numerical integration willinvariably be required. In Section 4, we analyze the limit of vanishing relaxation where analytical results areaccessible.It is worth noting at this point that the conductivity in Eq. 19 is not symmetrized. It follows from thedefinition in Eq. 6 that only the portion of the nonlinear conductivity that respects intrinsic permutationsymmetry is physical [3]. When permutations of Eq. 19 are properly accounted for, there will not only be anone-photon contribution associated to the resonance (cid:126) ω = ∆ (cid:15) cv but also (cid:126) ω = ∆ (cid:15) cv and so on. Likewisefor contributions associated with higher numbers of photons. Having this in mind, the formulas presentedfor the nonlinear conductivity in this work will nonetheless be left, for the most part, unsymmmetrized.Symmetrizing the conductivity is a trivial, if cumbersome, procedure and adds little to the discussion here.The terms described so far give a complete description for any two-band model that is used to describean insulator or a cold semiconductor. But for systems with charge carriers on the Fermi surface prior tooptical excitation, there is an additional contribution , σ αF (¯ ω , · · · , ¯ ω n ) = (cid:88) X = A,B,... (cid:88) p C pX (¯ ω , · · · , ¯ ω n ) F p ( α ) X (24)where the integrals have a different structure than before, F βα ··· α n X = (cid:90) d d k (2 π ) d (cid:88) a g α ...α n X ( A , (cid:15) ) aa ∂ β f a (25)It is evident from the presence of derivatives of Fermi functions in Eq. 25 that these integrals are deter-mined by the properties of the Fermi surface. More surprising is the absence of any frequency dependence.All dispersion in σ F comes from the C coefficients, which are the same for every two-band system. This leadsus to an important result: the dispersion of the contributions from the Fermi surface , those that dominatethe optical response of metals at low frequencies, is given by an universal family of functions of frequency (Eq. 24), obtained through linear combinations of the C X ’s. The particular linear combination observedis set by the integrals F X . Being dictated by Fermi surface properties, they are particularly dependent oncarrier concentration and can therefore be tuned by doping.In the following sections, Eqs. 20-25 are made explicit for the linear, second and third order conductivi-ties. Ready-to-use formulas are presented that reduce the calculation of the nonlinear conductivities to theevaluation of a minimal number of integrals over the FBZ (Eqs. 29-30, 37-39, 45-51). In linear order, the resonance-based decomposition of Eq. 19 falls into the familiar intra- and interbandseparation, σ βα (¯ ω ) = σ βα F (¯ ω ) + σ βα (¯ ω ) (26)with The validity of Eq. 24 has been checked by the authors only up to third order. βα F (¯ ω ) = − ie (cid:126) (cid:126) ¯ ω F βα A (27) σ βα (¯ ω ) = ie (cid:126) (cid:126) ¯ ω Π βα (¯ ω ) (28)The integrals required to obtain the linear response are well-known and particularly simple. There is oneintegral associated with an interband resonance and one related to the Fermi surface, F βα A ≡ (cid:90) d d k (2 π ) d (cid:88) a ∂ α (cid:15) a ∂ β f a (29)Π βα (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ab } (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (30)where the real and imaginary parts of products of Hermitian operators are symmetric and antisymmetric,respectively, under exchange of the band labels a and b ,Re { A ba B ab } = 12 ( A ba B ab + A ab B ba ) (31) i Im { A ba B ab } = 12 ( A ba B ab − A ab B ba ) (32) In second order, Eq. 19 gives σ βα α (¯ ω , ¯ ω ) = σ βα α F (¯ ω , ¯ ω ) + σ βα α (¯ ω , ¯ ω ) + σ βα α (¯ ω , ¯ ω ) (33)We now inspect each contribution separately. In the following, we make explicit the C coefficients inEq. 23 and the g functions in Eqs. 21 and 25 for n = 2.Starting with the Fermi surface contributions, σ βα α F (¯ ω , ¯ ω ) = ie (cid:126) (cid:126) ¯ ω F α α βB (34)The one-photon contributions, (cid:126) i e σ βα α (¯ ω , ¯ ω ) = 1 (cid:126) ¯ ω + (cid:126) ¯ ω Π α α β (¯ ω ) + (cid:126) ¯ ω + (cid:126) ¯ ω ( (cid:126) ¯ ω ) Π βα α (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω Π α βα (¯ ω ) (35)Finally, the two-photon contributions, (cid:126) i e σ βα α (¯ ω , ¯ ω ) = − (cid:126) ¯ ω + (cid:126) ¯ ω ( (cid:126) ¯ ω ) Π βα α (¯ ω + ¯ ω ) + (cid:126) ¯ ω + (cid:126) ¯ ω (cid:126) ¯ ω Π βα α (¯ ω + ¯ ω ) (36)In second order, there are two integrals from interband resonances and one from the Fermi surface, F βα α B ≡ (cid:90) d d k (2 π ) d (cid:88) a F α α a ∂ β f a (37)Π βα α (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Im {A βba A α ab } ( ∂ α ∆ (cid:15) ab ) (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (38)Π βα α (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Im {A βba A α ab ; α } (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (39)8here we made use of the generalized derivative notation introduced by Aversa and Sipe [17]: O ab ; α ≡ ( δ ab ∂ α − i ( A αaa − A αbb )) O ab .We chose, in Eq. 37, to label the integral over the Fermi surface with B instead of A to reserve thenotation F A to the purely intraband contribution. That is, F A is to be identified as the contribution thatwould be present even in a system where all interband transitions are neglected. It is however absent in thesecond order conductivity, due to time-reversal symmetry. It can easily be shown that the same is true forall even orders of the nonlinear conductivity.While the purely intraband contributions can be derived to any order without the need for a perturbationtheory as complex as the one discussed here, the inclusion of interband transitions brings new Fermi surfacecontributions when we move past linear order. In F βα α B in particular, we see the appearance of the quantity F βα a ≡ ∂ β A α aa − ∂ α A βaa which we recognize as the Berry curvature [27] and can be expressed as productsof the non-abelian Berry connection (Appendix A). It follows that the σ F contribution at second order is aprobe of the Berry curvature around the Fermi surface.Another curious result that follows from Eq. 37 is that we can expect the Fermi surface contributions tobe absent when the current is measured in the direction of the optical fields: σ βα α F = 0 when β = α = α .As a consequence, no Drude peaks should be found in the diagonal tensor elements of the second orderconductivity. Finally, we note that the contributions that probe the FBZ beyond the Fermi surface, the one- and two-photon contributions, are determined by the integrals in Eqs. 38 and 39 which are known in the literatureof the nonlinear optics of solids for their relation to the injection and shift currents in semiconductors [19],respectively.
A considerable jump in complexity occurs when we move to third order. Once again, we start with Eq. 19, σ βα α α (¯ ω , ¯ ω , ¯ ω ) = σ βα α α F (¯ ω , ¯ ω , ¯ ω )+ σ βα α α (¯ ω , ¯ ω , ¯ ω )+ σ βα α α (¯ ω , ¯ ω , ¯ ω )+ σ βα α α (¯ ω , ¯ ω , ¯ ω )(40)and write out all contributions explicitly.The Fermi surface contributions, σ βα α α F (¯ ω , ¯ ω , ¯ ω ) = ie (cid:126) (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω F βα α α A (41)The one-photon contributions, (cid:126) i e σ βα α α (¯ ω , ¯ ω , ¯ ω ) = − (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π α α α β (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω ( (cid:126) ¯ ω ) Π α α βα (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω Π α α α β (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π βα α α (¯ ω ) − (cid:126) ¯ ω ) (cid:126) ¯ ω Π α α α β (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω ) ( (cid:126) ¯ ω ) Π βα α α (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π α α α β (¯ ω ) + 2 (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π α α α β (¯ ω ) − ( (cid:126) ¯ ω ) (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π α α α β (¯ ω ) (42)with the abbreviations ¯ ω ij ≡ ¯ ω i + ¯ ω j and ¯ ω ≡ ¯ ω + ¯ ω + ¯ ω .The two-photon contributions, 9 i e σ βα α α (¯ ω , ¯ ω , ¯ ω ) =+ (cid:126) ¯ ω ( (cid:126) ¯ ω ) (cid:126) ¯ ω Π α α βα (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π α α α β (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω ( (cid:126) ¯ ω ) Π βα α α (¯ ω )+ 1( (cid:126) ¯ ω ) (cid:126) ¯ ω Π α α α β (¯ ω ) + (cid:126) ¯ ω ( (cid:126) ¯ ω ) ( (cid:126) ¯ ω ) Π βα α α (¯ ω ) − ( (cid:126) ¯ ω ) (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π α α α β (¯ ω ) (43)Finally, the three-photon contributions, (cid:126) i e σ βα α α (¯ ω , ¯ ω , ¯ ω ) = + (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω ( (cid:126) ¯ ω ) Π βα α α (¯ ω ) + (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π βα α α (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω ) ( (cid:126) ¯ ω ) Π βα α α (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π βα α α (¯ ω ) − (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω Π α α α β (¯ ω ) (44)In third order, there are six integrals from interband resonances and one from the Fermi surface, F βα α α A ≡ (cid:90) d d k (2 π ) d (cid:88) a ∂ α ∂ α ∂ α (cid:15) a ∂ β f a (45)Π βα α α (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ab ; α } ( ∂ α ∆ (cid:15) ab ) (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (46)Π βα α α (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ab } ( ∂ α ∂ α ∆ (cid:15) ab ) (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (47)Π βα α α (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ab } ( ∂ α ∆ (cid:15) ab ) ( ∂ α ∆ (cid:15) ab ) (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (48)Π βα α α (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ab ; α α } (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (49)Π βα α α (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ba A α ab A α ab } (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (50)Π βα α α (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba ; α A α ab ; α } (cid:126) ¯ ω − ∆ (cid:15) ab ∆ f ba (51)In these equations, all products of off-diagonal matrix elements of A appear in conjugate pairs (one A cv and one A vc ). This is required by gauge invariance: the nonlinear conductivity must be invariant under achange of phases of the Bloch functions ψ k a → e iθ k a ψ k a . For the same reason, generalized instead of regularderivatives of A must be used.The numerators of the Π αj integrals involve symmetric (Eq. 31) or antisymmetric (Eq. 32) combinationsof a product. Symmetric combinations occur for even orders of perturbation theory and antisymmetricones for odd orders. This structure is obtained by combining the contributions to the integrands from k and − k , then invoking time-reversal symmetry (see Appendix C for an example). This is an insightful andsometimes useful way to write the integrals. For instance, in evaluating diagonal tensor elements of thesecond order conductivity it is clear that the integral Π need not be considered, since it vanishes identically:Im {A βba A α ab } = 0 for β = α .The type of resonance-based analysis presented here could be continued for successively higher order non-linear conductivities. The number of integrals would continuously increase as more derivatives or generalizedderivatives are applied in the perturbation theory and more ways exist to arrange them in the numeratorsof the Π βα ...α n j integrals. 10 Relaxation-free limit
In Section 2.2, complex frequencies were introduced as a means to allow the two-band system to relax. Inthis section, the limit of vanishing relaxation and real frequencies is considered.As mentioned before, there are some subtleties in implementing this limit. For certain regions of fre-quency space, the nonlinear conductivity will diverge in the relaxation-free limit. Specifically, the nonlinearconductivity may diverge in the relaxation-free limit when a subset of the optical frequencies adds to zero.This includes the case of DC fields ( ω i = 0) or DC currents ( ω + · · · + ω n = 0).Fundamentally, this is due to the existence of regions of the FBZ where more than one of the resonanceconditions is obeyed simultaneously (two or more denominators in the symmetrized version of Eq. 18 areresonant at those k -points). Geometrically, this can be visualized by representing in the FBZ the set ofpoints where the conditions ω − ∆ (cid:15) k cv = 0 or ω + · · · + ω i − ∆ (cid:15) k cv = 0 ( i = 2 , ..., n ) and their permutationsare verified and checking for overlapping regions (for example, in Fig. 2b for there to be the possibility of adivergence two of the contours would have to cross or altogether overlap).When using the equations of Section 3, this difficulty is realized in the relaxation-free limit of Eq. 23,lim γ → + σ αi (¯ ω , · · · , ¯ ω n ) = lim γ → + (cid:88) j, p C pij (¯ ω , · · · , ¯ ω n ) Π p ( α ) j (¯ ω + · · · +¯ ω i ) = (cid:88) j, p C pij ( ω , · · · , ω n ) lim γ → + Π p ( α ) j (¯ ω + · · · +¯ ω i )(52)The validity of this passage hangs on the convergence of the C coefficients in the relaxation-free limit,which in turn relies on there existing no cancelling subset of frequencies in { ω , . . . , ω n } .Henceforth, we will assume this to be the case. We emphasize that the validity of the equations in Section 3always holds, even when in these regions of frequency space. However, the analysis of the relaxation-freelimit in these situations, which include several nonlinear optical effects of interest (photovoltaic effects, DCfield-induced second order response, electro-optic effects,...), albeit possible, requires greater care and willbe left to future communications. Indeed, the study of the divergences of the nonlinear optical conductivityis a topic of current research interest [28–30]. From Eq. 52, the evaluation of the relaxation-free limit of the nonlinear conductivity has been reduced toevaluating the relaxation-free limit of the integrals Π αj . This is a direct application of Eq. 17,lim γ → + Π αj (¯ ω ) = π (cid:0) H αj ( ω ) − i I αj ( ω ) (cid:1) (53)where I αj ( ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b g αj ( A , (cid:15) ) ab ∆ f ba δ ( (cid:126) ω − ∆ (cid:15) ab ) (54) H αj ( ω ) ≡ − (cid:90) d d k (2 π ) d (cid:88) a,b g αj ( A , (cid:15) ) ab (cid:126) ω − ∆ (cid:15) ab ∆ f ba (55)As an example, let us consider the linear conductivity. Ignoring Fermi surface contributions,lim γ → + σ βα (¯ ω ) = ie ω (cid:16) H βα ( ω ) − i I βα ( ω ) (cid:17) (56)where I βα ( ω ) = (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ab } ∆ f ba δ ( (cid:126) ω − ∆ (cid:15) ab ) (57) H βα ( ω ) = − (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ab } (cid:126) ω − ∆ (cid:15) ab ∆ f ba (58)11he integral H can be computed through Eq. 58 or from applying a Hilbert transform to I . This istrue in general, as it is easily proven from Eqs. 54 and 55 that H αj ( ω ) = − π − (cid:90) + ∞−∞ I αj ( ω (cid:48) ) ω (cid:48) − ω dω (cid:48) (59)The entire linear response of an insulator or a cold semiconductor is obtained, for negligible relaxation,from a calculation of the function in Eq. 57 and its Hilbert transform. This is by no means new informationin linear optics: the integral in Eq. 57 is nothing more than Fermi’s golden rule and the Hilbert transformfrom I to H equates to the usual derivation of refraction from absorption by the Kramers-Kr¨onig relations.Here, however, this method has been generalized to any order.The real part of the nonlinear conductivity is given by the I integrals and the imaginary part by the H integrals . This can be seen by noting that the coefficients C that multiply the integrals in Eq. 52 are purelyimaginary in the relaxation-free limit, while I ( ω ) and H ( ω ) are always real.The association of the real part of the nonlinear conductivity with the I integrals is to be expected if onerecalls that the real part is the one responsible for optical absorption. If the optical frequencies (and theirharmonics) of the light incident in a cold semiconductor all lie below the gap, there must be no absorptionand therefore the real part of the nonlinear conductivity is required to vanish. This is guaranteed to be thecase for the integrals in Eq. 54.We conclude from the analysis of this section that, if we wish to evaluate the real part of a nonlinearconductivity, we have only to evaluate simple integrals with the form of Eq. 54 for the appropriate g functions,listed in Eqs. 29-30, 37-39, 45-51. These contributions are well localized in the FBZ and run over regionsthat satisfy the resonance conditions.The imaginary part of a nonlinear conductivity follows directly by performing Hilbert transforms of theintegrals computed for the real part. We note that we are not performing a Hilbert transform of the entirenonlinear conductivity as it is traditionally done in the nonlinear Kramers-Kr¨onig relations [31]. We arenot even applying a Hilbert transform to the different contributions in Eq. 19 to move from their real toimaginary parts . The passage from the real to imaginary parts of the nonlinear optical conductivity is moresubtle and is made directly through the integrals in Eq. 53. In this way, the imaginary part, which unlikethe real part is not described by any restricted region of the FBZ, is made more accessible than what wouldperhaps be expected.Finally, we point out that for Fermi surface contributions there is no need for the relaxation-free limit,since the F X integrals are frequency independent. To attempt an analytical calculation it is necessary onlyto set the temperature to zero, in which case F βα ··· α n X = − (cid:90) d d k (2 π ) d (cid:88) a g α ...α n X ( A , (cid:15) ) aa ∂ β (cid:15) a δ ( µ − (cid:15) a ) (60)This integral has a similar structure to the integrals in Eq. 54, with the chemical potential taking therole of a frequency.For the previous integrals (Eq. 54), it is equally helpful to consider the system to be at zero temperatureto make calculations tractable. However, once the nonlinear conductivity is calculated for T = 0, it ispossible to obtain the answers at finite temperature with a few tricks [23]. It is possible to quickly arrive at finite temperature results from zero temperature calculations of the nonlinearconductivity by making use of a relation first introduced in [23]. The key idea is to recognize that the entirechemical potential dependence of the Π αi (¯ ω ) and F αX integrals stems from the Fermi-Dirac distribution.Starting with the former type of integrals, we make all arguments of the distribution explicit, f ( (cid:15), µ, T ) ≡ e ( (cid:15) − µ ) /k B T + 1 (61) Once again, excluding the regions of frequency space where subsets of frequencies add to zero. This would not work as the real and imaginary parts of each individual contribution are not simply related by a Hilberttransform.
12t zero temperature, f ( (cid:15), µ,
0) = Θ( µ − (cid:15) ) (62)The finite temperature and zero temperature distributions are related by f ( (cid:15), µ, T ) = (cid:90) + ∞−∞ f ( µ (cid:48) , µ, T ) δ ( µ (cid:48) − (cid:15) ) dµ (cid:48) = (cid:90) + ∞−∞ f ( µ (cid:48) , µ, T ) ∂ µ (cid:48) f ( (cid:15), µ (cid:48) , dµ (cid:48) (63)Since no other objects in Π αi (¯ ω ) depend on the chemical potential, Eq. 63 translates directly to σ βα ...α n i (¯ ω , . . . , ¯ ω n ; µ, T ) = (cid:90) + ∞−∞ f ( µ (cid:48) , µ, T ) ∂ µ (cid:48) σ βα ...α n i (¯ ω , . . . , ¯ ω n ; µ (cid:48) , dµ (cid:48) (64)where we made the chemical potential and temperature dependencies explicit.By a similar argument, σ F is proven to also obey Eq. 64. Joining all contributions, we arrive at a generalrelation between the zero and finite temperature nonlinear conductivity, σ βα ...α n (¯ ω , . . . , ¯ ω n ; µ, T ) = (cid:90) + ∞−∞ f ( µ (cid:48) , µ, T ) ∂ µ (cid:48) σ βα ...α n (¯ ω , . . . , ¯ ω n ; µ (cid:48) , dµ (cid:48) = 1 T (cid:90) + ∞−∞ (1 − f ( µ (cid:48) , µ, T )) f ( µ (cid:48) , µ, T ) σ βα ...α n (¯ ω , . . . , ¯ ω n ; µ (cid:48) , dµ (cid:48) (65)by partial integration, using σ ( µ → −∞ ) = 0.The effect of having a finite temperature is to probe the nonlinear conductivity at different values of thechemical potential. Eq. 65 is a kind of weighted average given by a distribution that is centered and peakedat the Fermi level and has a width of the order k B T .We conclude by making some observations on systems with electron-hole symmetry for which the chemicalpotential dependence of σ αi can be specially simple.The real part of the one-photon contribution is given by the I αj ( ω ) integrals. At T = 0 and assumingelectron-hole symmetry, their integrand is proportional to,∆ f vc = Θ( µ − (cid:15) v ) − Θ( µ − (cid:15) c ) = Θ(∆ (cid:15) cv − | µ | ) → Θ( (cid:126) | ω | − | µ | ) (66)where the last step uses the Dirac delta in Eq. 54. The final Heaviside step function is independent of k andcan be pulled out of the I αj ( ω ) integrals: we have therefore an universal dependence on chemical potential,with frequency dependent coefficients set by the material, the inverse situation to what happens in σ F . Forgapped systems, we could more generally write Θ( (cid:126) | ω | − Max(2 | µ | , ∆)), where Max(2 | µ | , ∆) is the “effectivegap”.Similarly, for the i -photon contribution,∆ f vc → Θ( (cid:126) | ω + · · · + ω i | − Max(2 | µ | , ∆)) (67)This is exemplified by the expression in Eq. 1. The real part of the nonlinear conductivity of electron-hole symmetric systems is always given in termsof step functions in the chemical potential , aside from the low frequency behaviour where Fermi surfacecontributions dominate. These step-like variations are due to Pauli blocking. Depending on the dispersionof the coefficients of the step functions, the steps can be found to vary in magnitude and even be absent orinfinite.On the other hand, the imaginary part is obtained by applying a Hilbert transform to I αi and differentsystems will give different functions of the chemical potential. No universal behaviour can be predicted, inthis context, for the imaginary part.The corresponding finite temperature result follows from Eqs. 66 and 64, (cid:90) + ∞−∞ f ( µ (cid:48) , µ, T ) ∂ µ (cid:48) Θ( (cid:126) | ω | − | µ | ) dµ (cid:48) = f ( − (cid:126) | ω | / , µ, T ) − f ( (cid:126) | ω | / , µ, T ) (68)Introducing a finite temperature smooths the step functions observed when varying the chemical potential.13 A simplistic model for direct gap semiconductors
A detailed framework for computations of the nonlinear conductivity has been described and it is now leftto demonstrate it with some simple models. We consider first the case of direct gap semiconductors (3D). Ofcourse, these materials can be not only of great interest, but of great complexity. If we focus on frequencieswhose resonances lie close to the gap, however, the band structure may be approximated by a single pair ofisotropic parabolic bands with effective masses, (cid:15) k c = (cid:126) k m c + ∆2 (69) (cid:15) k v = − (cid:126) k m v − ∆2 (70)with k = (cid:113) k x + k y + k z .The trouble now lies in describing the non-abelian Berry connection, whose derivation relies on a knowl-edge of the energy eigenstates. A common approach that was adopted in early attempts at describing thenonlinear properties of solids [3, 20, 32, 33], and is still in use in the literature [34], is to consider nearlyconstant dipole matrix elements. In the language of this work, this means assuming not only that the bandsare parabolic near the gap but that the gauge field A β is constant over the region of the FBZ that is beingprobed by the optical fields. This overly simplistic description of semiconductors has at least the merit ofproviding quick derivations and a first intuition on the dispersion of the nonlinear conductivity near the gapof a semiconductor.Somewhat paradoxically, A β is sometimes taken to be isotropic as well [20, 32, 33]. This is clearly notpossible since for a vector field to be isotropic and constant it either vanishes or it is not uniquely defined overthe region of k -space that is of interest. Since no degeneracies are being considered, the non-abelian Berryconnection A β must be well-defined. We must then break the isotropy and consider a particular directionthat the field is aligned with, say along the x axis, A βcv = A δ βx (71)Eqs. 69, 70 and 71 define the model for which we shall compute optical conductivities. For simplicity, weconsider the incident optical field to be linearly polarized along x and derive the response along the samedirection.For the undoped system, the linear conductivity follows from the evaluation of the integral in Eq. 57, I xx ( ω ) = |A| Sign( ω ) ρ JDOS ( (cid:126) | ω | ) (72)where JDOS stands for joint density of states, ρ JDOS ( (cid:126) | ω | ) = (cid:90) d k (2 π ) δ ( (cid:126) | ω | − ∆ (cid:15) k cv ) = 1 √ π (cid:126) m / r (cid:112) (cid:126) | ω | − ∆ Θ ( (cid:126) | ω | − ∆) (73)and m − r ≡ m − c + m − v .From Eq. 28,Re ( σ xx ( ω )) = π e (cid:126) (cid:126) ω I xx ( ω ) = e √ π (cid:126) m / r |A| (cid:126) | ω | (cid:112) (cid:126) | ω | − ∆ Θ ( (cid:126) | ω | − ∆) (74)The real part of the linear conductivity is depicted in Fig. 1a. The dispersion observed near the gap isdirectly related to the well-known square-root dependence of optical absorption on frequency [3].Since the joint density of states vanishes at the gap, there are no steps in the dispersion of the linearconductivity until the system is doped and an effective gap is set at higher energies (Fig. 1a). The linearconductivity is then given by The Fermi surface contribution gives the Drude conductivity, whose real part consists of a Dirac delta at zero frequency inthe relaxation-free limit. Since we excluded zero frequencies from our analysis (see the discussion in Section 4) the real part ofthe Fermi surface contribution is neglected in Eqs. 75, 90 and 93. = μ = .075 Ry0.00 0.05 0.10 0.15 0.20 0.25024681012 ℏω ( Ry ) R e [ σ xx ( ω ) ] ( S m - ) (a) μ = μ = .075 Ry0.00 0.05 0.10 0.15 0.20 - - ℏω ( Ry ) R e [ σ xxxx ( ω , ω , ω ) ] ( - S m V - ) (b) Figure 1:
Real part of the (a) linear and (b) third order conductivity for a system of parabolic bands. The energiesare given in Rydberg’s (1 Ry (cid:39) . eV ). The chosen parameters are: m c = m v = 0 . m e , A = 6 . . .
15 Ry in the dopedcase. The value chosen for A is derived from the value of E p in [33]. Re ( σ xx ( ω )) = Re ( σ xx ( ω )) = e √ π (cid:126) m / r |A| (cid:126) | ω | (cid:112) (cid:126) | ω | − ∆ Θ ( (cid:126) | ω | − ∆ eff ) (75)with the effective gap, ∆ eff = (cid:40) m c m r ( µ − ∆ /
2) + ∆ µ > m v m r ( − µ − ∆ /
2) + ∆ µ ≤ k dependence of A is neglected, the generalized derivatives are taken to be identically zero and half of the integrals to beevaluated at third order vanish, I xxxx ( ω ) = I xxxx ( ω ) = I xxxx ( ω ) = 0 (77)The remaining integrals are straightforwardly computed (Eqs. 47, 48 and 50), I xxxx ( ω ) = (cid:126) m r |A| ρ JDOS ( (cid:126) | ω | ) = 1 √ π (cid:126) m / r |A| (cid:112) (cid:126) | ω | − ∆ Θ ( (cid:126) | ω | − ∆ eff ) (78) I xxxx ( ω ) = √ π (cid:126) m / r |A| Sign( ω ) ( (cid:126) | ω | − ∆) / Θ ( (cid:126) | ω | − ∆ eff ) (79) I xxxx ( ω ) = |A| Sign( ω ) ρ JDOS ( (cid:126) | ω | ) = 1 √ π (cid:126) m / r |A| Sign( ω ) (cid:112) (cid:126) | ω | − ∆ Θ ( (cid:126) | ω | − ∆ eff ) (80)Direct substitution of these integrals in the relaxation-free limit of Eqs. 40-44 immediately gives a generalexpression for the third order conductivity. For brevity, we make the expression explicit only for the specialcase of harmonic generation ( ω = ω = ω = ω ), σ xxxx ( ω, ω, ω ) = σ xxxx ( ω, ω, ω ) + σ xxxx ( ω, ω, ω ) + σ xxxx ( ω, ω, ω ) (81)15 F (a) σ σ F σ σ (b) Figure 2: (a) Cross section of the Dirac cone [35] with the dark lines showing the linear energy-momentum relationin Eqs. 85 and 86. The arrows stand for incident photons, all of which are assumed to have the same energy (thecolor scheme serves only to differentiate the different contributions) but probe different regions of the FBZ dependingon the number of photons involved. (b) Contours in the FBZ where the resonance conditions that define the differentcontributions in Eq. 19 are met. For the two-dimensional crystal of monolayer graphene, they consist of circlescentered at the Dirac point. The red circle is the Fermi surface. where the one-, two- and three-photon contributions are σ xxxx ( ω, ω, ω ) = πie (cid:126) (cid:18) − I xxxx ( ω )4 ( (cid:126) ω ) − I xxxx ( ω )6 ( (cid:126) ω ) + 3 I xxxx ( ω )4 (cid:126) ω (cid:19) (82) σ xxxx ( ω, ω, ω ) = πie (cid:126) I xxxx (2 ω )3 ( (cid:126) ω ) (83) σ xxxx ( ω, ω, ω ) = πie (cid:126) (cid:18) I xxxx (3 ω )4 ( (cid:126) ω ) − I xxxx (3 ω )2 ( (cid:126) ω ) − I xxxx (3 ω )4 (cid:126) ω (cid:19) (84)The Fermi surface contributions are absent at third order for a quadratic dispersion relation ( ∂ (cid:15) = 0 inEq. 45).The real part of the third order conductivity is represented in Fig. 1b. Three features are clearly seen forincreasing frequency when the three-, two- and one-photon contributions, respectively, become relevant. Forthe doped case, the usual steps are observed, stemming from the Heaviside theta functions in Eqs. 78, 79and 80, which in turn are manifestations of Pauli blocking of transitions below the Fermi level. In theundoped case, analogous to linear order, the steps are suppressed by the vanishing joint density of states,giving rise to pronounced, but less abrupt features.Once again, the imaginary part is not captured by the current description. As a final example, we turn to monolayer graphene [35, 36] whose nonlinear optical properties have been asubject of active research since the pioneering works of Mikhailov et al. [37–39]. A widely adopted two-bandtight binding model provides the band structure and non-abelian Berry connection near the Dirac point [22], (cid:15) k c = + v F (cid:126) k (85) (cid:15) k v = − v F (cid:126) k (86) A βcv = ( ˆz × k ) β k (87)with ˆz = (0 , , k = (cid:113) k x + k y and v F (cid:39) m/s as the Fermi velocity. The energy-momentum relation isdepicted in Fig. 2a together with the contours defined by photon resonances in Fig. 2b16 = μ = - - ℏω / t R e [ σ xxxx ( ω , ω , ω ) ] ( - S m V - ) (a) σ σ σ - - - - ℏω / t R e [ σ xxxx ( ω , ω , ω ) ] ( - S m V - ) (b) Figure 3:
Real part of the third order conductivity of monolayer graphene as a function of the optical frequency,normalized to the tight binding parameter t (cid:39) eV [35], near the Dirac point. In (a), both undoped and dopedgraphene are shown. In (b), the different contributions to the conductivity (see Eq. 19 and Fig. 2b) are representedfor µ = 0 . t . In the inset, it is possible to discern a small feature in the one-photon contribution. We leave the integral evaluation to an appendix (Appendix D) and proceed to inspect the analyticexpressions for the conductivities. In linear order, symmetry reduces the number of independent tensorelements to one, σ xy ( ω ) = σ yx ( ω ) = 0 (88) σ yy ( ω ) = σ xx ( ω ) = σ xxF ( ω ) + σ xx ( ω ) (89)with σ xxF ( ω ) = 4 i σ π µ (cid:126) ω (90) σ xx ( ω ) = σ (cid:18) Θ ( (cid:126) | ω | − | µ | ) + iπ ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (91)For | ω | (cid:29) | µ | , we find the universal conductivity of graphene σ xx ( ω ) = σ = e / (cid:126) [40].Fortunately, in this model the Hilbert transforms converge (Eqs. 154 and 174 of Appendix D) andexpressions for the imaginary part are obtainable, as included in Eq. 91. The case of graphene showcasesa general characteristic of the optical response that stems from the Hilbert transforms relating I and H : discontinuities in the real part (as those introduced by Heaviside theta functions) translate into logarithmicdivergences in the imaginary part and vice-versa.Monolayer graphene is inversion symmetric so it has no second order response.Evaluation of the integrals in Eqs. 45-51 gives a third order conductivity (Appendix D) in agreementwith the results in the literature [22]. The third order response of graphene has already been the subject ofseveral theoretical studies [22, 23, 37, 38, 41] and the form of the third order conductivity in the independentelectron approximation is known [22, 42]. Here, we merely aim to exemplify how the equations of Section 3provide a straightforward derivation and how the structure of the nonlinear conductivity emerges naturallyin this framework.For these purposes, it is sufficient to consider third harmonic generation where symmetry constraintsreduce the tensor to a single independent component σ xxyy ( ω, ω, ω ) = σ xyxy ( ω, ω, ω ) = σ xyyx ( ω, ω, ω ) = σ xxxx ( ω, ω, ω ) / σ xxxx ( ω, ω, ω ) = σ xxxxF ( ω, ω, ω ) + σ xxxx ( ω, ω, ω ) + σ xxxx ( ω, ω, ω ) + σ xxxx ( ω, ω, ω ) (92)17 = μ = - - ℏω / t I m [ σ xxxx ( ω , ω , ω ) ] ( - S m V - ) (a) σ F σ σ σ - - - - ℏω / t I m [ σ xxxx ( ω , ω , ω ) ] ( - S m V - ) (b) Figure 4:
Imaginary part of the third order conductivity of monolayer graphene as a function of the optical frequency,normalized to the tight binding parameter t (cid:39) eV [35], near the Dirac point. In (a), both undoped and dopedgraphene are shown. In (b), the different contributions to the conductivity (see Eq. 19 and Fig. 2b) are representedfor µ = 0 . t . In the inset, it is possible to discern a small feature in the one-photon contribution. with σ xxxxF ( ω, ω, ω ) = 24 i (cid:126) C π | µ | ω (93) σ xxxx ( ω, ω, ω ) = − C ω (cid:18) Θ ( (cid:126) | ω | − | µ | ) + iπ ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) + 71 i (cid:126) C π | µ | ω (94) σ xxxx ( ω, ω, ω ) = + 64 C ω (cid:18) Θ (2 (cid:126) | ω | − | µ | ) + iπ ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) − i (cid:126) C π | µ | ω (95) σ xxxx ( ω, ω, ω ) = − C ω (cid:18) Θ (3 (cid:126) | ω | − | µ | ) + iπ ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) − i (cid:126) C π | µ | ω (96)The result is plotted in Figs. 3 and 4. Inspection of Eqs 94-96 and the resulting curves shows thatthe different contributions are, as expected, dominant in the specific regions of the spectrum where therelevant photon frequencies match the gap. Once again the real part of the conductivity is marked by jumpdiscontinuities and the imaginary part by the associated logarithmic divergences. The real parts of the one-,two- and three-photon contributions are zero below the gap, while the imaginary parts tend to zero in theDC limit. In this limit, the Fermi surface contribution (Eq. 93) is dominant.A curious special case is that of undoped graphene. A hasty look through the expressions for the thirdorder conductivity might lead one to think that it diverges when the chemical potential is set to zero. Butwhen summing all the contributions (Eqs. 93-96) to the third order conductivity the different factors of µ − cancel (Eq. 1).Nonetheless, there have been some concerns with this divergence in the literature. In the work of Chenget al. [23] a divergence of such nature was found. But, as pointed out by the authors, it appears only whenintroducing two distinct damping parameters, one for intraband and one for interband resonances. Thissuggests that such a divergence is an artifact of the phenomenological description of relaxation that wasadopted in their work. In the method used here, there is a single parameter γ describing relaxation, but thisparameter might be made different at low and high frequencies (as pointed out before, any even function γ = γ ( ω ) will do) and still no divergences for µ = 0 are found.Similarly, for the Boltzmann equation approach adopted in the pioneering works of [37, 38], the thirdorder conductivity was found to diverge when µ →
0. The reason is transparent in our framework: theBoltzmann equation approach takes only the contribution from the Fermi surface into account (purelyintraband transitions) and this one can be seen from Eq. 93 to indeed diverge for a zero chemical potential.18oltzmann equation methods are appropriate to describe the low frequency nonlinear conductivity as longas there is a gap, either a real or an effective one set by Pauli blocking, below which all contributions otherthan the ones defined by the Fermi surface become increasingly negligible. For undoped graphene however,this is never the case and the remaining contributions in Eq. 92 must be considered.Unfortunately, a real difficulty can still be found for undoped systems with zero gap in the DC limit.Even though the µ = 0 third order conductivity of graphene is well defined for any finite frequency, it showsa frequency dependence ω − (see Eq. 2 of [22]) diverging when ω →
0. A similar point could be made aboutthe model in the previous section if we close the gap (∆ eff = ∆ = 0 in Eqs. 78-80). This result holds evenfor a tight binding model defined over the entire FBZ and is presumably due to a breakdown in perturbationtheory.
A decomposition of the nonlinear conductivity based on the possible resonances between optical frequenciesand band energies is useful, both from the perspective of practical calculations as well as in developing anphysical intuition and an understanding of the structure of this response function. This was exemplified withan extensive discussion of the nonlinear optical response of the time-reversal symmetric two-band crystal.Calculations of the first, second and third order conductivity have been shown to reduce to the evaluationof a small number of integrals over the FBZ. Evaluation of these integrals may be performed numerically butalso analytically if one considers the relaxation-free limit. In this limit, obtaining the real part of a nonlinearconductivity is no harder than performing a series of Fermi golden rule computations. The imaginary partfollows by applying Hilbert transforms directly to the results of these computations, without need of thenonlinear Kramers-Kr¨onig relations.A known difficulty in implementing the relaxation-free limit is that the nonlinear conductivity is thensingular for extended regions of the frequency plane. These singularities are behind many of the moreinteresting optical effects that have been studied in recent years, such has the shift, injection and jerkcurrents [29, 43]. For brevity, we did not address these effects here, but the equations in Section 3 providea natural starting point for the derivation of expressions for small but finite γ .The effect of temperature was considered and some simple models solved for the purpose of illustratingthe method. The third order conductivity of monolayer graphene was shown to be well-defined in the limit µ → Acknowledgments
The authors acknowledge funding from Funda¸c˜ao da Ciˆencia e Tecnologia (FCT), under the COMPETE2020 program in FEDER component (European Union), through projects POCI-01-0145-FEDER-028887,UID/FIS/04650/2019 and M-ERA-NET2/0002/2016. The work of G. B. V. and D. J. P. is supported bythe grants PD/BD/140891/2018 and PD/BD/135019/2017, respectively.19 ppendix A: Commuting position operators
It is implicitly assumed throughout this work that position operators commute. This may seem a trivialstatement, derived from the basics of quantum mechanics, but commutation relations tend to be brokenupon band truncation and when a finite band model is used to describe a system, defined perhaps by thespecification of a set of Bloch states ψ k a and bands (cid:15) k a , there is no reason to assume that the following willbe true, (cid:2) ˆ r α , ˆ r β (cid:3) = − (cid:104) ˆ D α , ˆ D β (cid:105) = 0 (97)Still, this particular commutation relation can be shown to hold on very general grounds, (cid:104) ˆ D α , ˆ D β (cid:105) k ab = − i ∂ α A β k ab + i ∂ β A α k ab − (cid:2) A α , A β (cid:3) k ab = ∂ α (cid:0)(cid:10) u k a (cid:12)(cid:12) ∂ β u k b (cid:11)(cid:1) − ∂ β ( (cid:104) u k a | ∂ α u k b (cid:105) ) + (cid:88) c (cid:0) (cid:104) u k a | ∂ α u k c (cid:105) (cid:10) u k c (cid:12)(cid:12) ∂ β u k b (cid:11) − (cid:10) u k a (cid:12)(cid:12) ∂ β u k c (cid:11) (cid:10) u k c (cid:12)(cid:12) ∂ β u k b (cid:11)(cid:1) = (cid:10) ∂ α u k a (cid:12)(cid:12) ∂ β u k b (cid:11) − (cid:10) ∂ β u k a (cid:12)(cid:12) ∂ α u k b (cid:11) − (cid:88) c (cid:0) (cid:104) ∂ α u k a | u k c (cid:105) (cid:10) u k c (cid:12)(cid:12) ∂ β u k b (cid:11) − (cid:10) ∂ β u k a (cid:12)(cid:12) u k c (cid:11) (cid:104) u k c | ∂ α u k b (cid:105) (cid:1) = 0 (98)with the last step making use of the closure relation, (cid:88) c | u k c (cid:105) (cid:104) u k c | = ˆ (99)Eq. 99 is the essential assumption in the construction of the finite-band model that ensures the validityof Eq. 97. With the commutation of position operators on firm grounds, some of its consequences can bediscussed.A curious and sometimes useful way to rewrite the statement of Eq. 97 is (cid:104) ˆ D α , ˆ A β (cid:105) ab = ∂ β A αab (100)where we dropped the k label.From Eq. 97, various other identities can be derived. Expanding the commutator, ∂ α A βab − ∂ β A αab = i (cid:2) A α , A β (cid:3) ab (101)In the particular case that a = b , it relates the abelian Berry curvature to the off-diagonal matrix elementsof the Berry connection, F αβ ≡ ∂ α A βaa − ∂ β A αaa = i (cid:88) c (cid:54) = a,b (cid:0) A αac A βca − A βac A αca (cid:1) (102)For a two-band model, Eq. 101 with a (cid:54) = b relates “generalized derivatives” of Berry connections, A βab ; α = A αab ; β (103)A more sophisticated identity on second order “generalized derivatives” can also be derived, A βab ; α α − A βab ; α α = 2 (cid:16) A α ba A βab A α ab − A α ba A βab A α ab (cid:17) (104)which translates into a relation between Π and Π which are no longer independent,Π βα α α (¯ ω ) − Π βα α α (¯ ω ) = 2 (cid:16) Π βα α α (¯ ω ) − Π βα α α (¯ ω ) (cid:17) (105)20 ppendix B: Integral identities In the derivation of the expressions in Sections 3.1-3.3, identities were used involving the parity and thederivatives of the integrals in the Eqs. 30, 38-39, 46-51. These identities are helpful in manipulating generalexpressions and are summarized here. Their use is exemplified in Appendix C.The first set of identities concerns the derivatives of some Π αi integrals. In second order, dd ¯ ω Π βα α (¯ ω ) = Π βα α (¯ ω ) + Π α βα ( − ¯ ω ) + Π βα α B (¯ ω ) (106)with Π βα α B (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Im {A βba A α ab } ( ∂ α ∆ f ba ) (cid:126) ¯ ω − ∆ (cid:15) ab (107)which vanishes in the absence of a Fermi surface.In third order, dd ¯ ω Π βα α α (¯ ω ) = Π βα α α (¯ ω ) + Π βα α α (¯ ω ) + Π βα α α B (¯ ω ) (108) dd ¯ ω Π βα α α (¯ ω ) = Π βα α α (¯ ω ) − Π α βα α ( − ¯ ω ) + Π βα α α (¯ ω ) + Π βα α α C (¯ ω )= Π βα α α (¯ ω ) − Π α βα α ( − ¯ ω ) + Π βα α α (¯ ω ) + Π βα α α C (¯ ω ) (109)with Π βα α α B (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ab ; α } ( ∂ α ∆ f ba ) (cid:126) ¯ ω − ∆ (cid:15) ab (110)Π βα α α C (¯ ω ) ≡ (cid:90) d d k (2 π ) d (cid:88) a,b Re {A βba A α ab } ( ∂ α ∆ (cid:15) ab ) ( ∂ α ∆ f ba ) (cid:126) ¯ ω − ∆ (cid:15) ab (111)which vanish in the absence of a Fermi surface.Additionally, from Eq. 109 it also follows,Π βα α α (¯ ω ) − Π βα α α (¯ ω ) − Π α βα α ( − ¯ ω ) + Π α βα α ( − ¯ ω ) = Π βα α α C (¯ ω ) − Π βα α α C (¯ ω ) (112)The second set of identities, unlike the first, is reliant on time-reversal symmetry and consists on simplystating the parity of the integrals in Eqs. 30, 38-39, 46-51,Π βα ( − ¯ ω ) = Π βα (¯ ω ) (113)Π βα α ( − ¯ ω ) = Π βα α (¯ ω ) (114)Π βα α ( − ¯ ω ) = − Π βα α (¯ ω ) (115)Π βα α α ( − ¯ ω ) = − Π βα α α (¯ ω ) (116)Π βα α α ( − ¯ ω ) = − Π βα α α (¯ ω ) (117)Π βα α α ( − ¯ ω ) = Π βα α α (¯ ω ) (118)Π βα α α ( − ¯ ω ) = Π βα α α (¯ ω ) (119)Π βα α α ( − ¯ ω ) = Π βα α α (¯ ω ) (120)Π βα α α ( − ¯ ω ) = Π βα α α (¯ ω ) (121)21inally, there is a third set of identities that relate different tensor elements of a given integral, F βα A = F α βA (122)Π βα (¯ ω ) = Π α β (¯ ω ) (123) F βα α B = − F α βα B (124)Π βα α (¯ ω ) = − Π α βα (¯ ω ) (125)Π βα α (¯ ω ) = Π βα α (¯ ω ) (126) F βα α α A = F α βα α A = F α α βα A = F α α α βA (127)Π βα α α (¯ ω ) = Π βα α α (¯ ω ) (128)Π βα α α (¯ ω ) = Π α βα α (¯ ω ) = Π βα α α (¯ ω ) (129)Π βα α α (¯ ω ) = Π α βα α (¯ ω ) = Π βα α α (¯ ω ) (130)Π βα α α (¯ ω ) = Π βα α α (¯ ω ) (131)Π βα α α (¯ ω ) = Π α βα α (¯ ω ) = Π βα α α (¯ ω ) = Π α α βα (¯ ω ) (132)Π βα α α (¯ ω ) = Π α βα α (¯ ω ) = Π βα α α (¯ ω ) = Π α α βα (¯ ω ) (133)These can be seen by inspection of the integrals, with Eqs. 126, 128, 131 and 133 as a consequence ofEq. 103.These symmetries reduce the number of independent tensor elements that are required for a nonlinearconductivity calculation. Appendix C: Derivation of the second order conductivity
As an illustration of the general scheme by which the expressions in Sections 3.1-3.3 were derived, the caseof the second order conductivity is treated here in detail. The aim is to arrive at the results of Section 3.2,starting with the general expression for a nonlinear conductivity. The derivation of third order conductivity,albeit considerably lengthier, follows along the same lines.We start with Eq. 18, taken for the case n = 2, σ βα α (¯ ω , ¯ ω ) = (cid:90) d d k (2 π ) d (cid:88) a,b J β k ba (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab (cid:20) r α , (cid:126) ¯ ω − ∆ (cid:15) k ◦ [ r α , ρ ] (cid:21) k ab (134)Substituting for the current matrix elements, Eq. 15, and using ˆ r α = i ˆ D α , σ βα α (¯ ω , ¯ ω ) = − e (cid:126) (cid:90) d d k (2 π ) d (cid:88) a,b (cid:104) ˆ D β , H (cid:105) k ba (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab (cid:20) r α , (cid:126) ¯ ω − ∆ (cid:15) k ◦ [ r α , ρ ] (cid:21) k ab = e (cid:126) (cid:90) d d k (2 π ) d (cid:88) a ∂ β (cid:15) k a (cid:126) ¯ ω + (cid:126) ¯ ω (cid:20) D α , (cid:126) ¯ ω − ∆ (cid:15) k ◦ [ D α , ρ ] (cid:21) k aa − i e (cid:126) (cid:90) d d k (2 π ) d (cid:88) a (cid:54) = b A β k ba ∆ (cid:15) k ab (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab (cid:20) D α , (cid:126) ¯ ω − ∆ (cid:15) k ◦ [ D α , ρ ] (cid:21) k ab (135)where diagonal and off-diagonal current matrix elements were separated.To expand the condensed expression in Eq. 135 into more explicit formulae, Eq. 14 is applied iteratively,starting with [ ˆ D α , ˆ ρ ] k ab = δ ab ∂ α f k a − i A α k ab ∆ f k ba (136)22here δ ab is the Kronecker delta.In the intermediate steps, the definition of the Hadamard product is used. For instance, (cid:18) (cid:126) ¯ ω − ∆ (cid:15) k ◦ [ D α , ρ ] (cid:19) k aa = [ D α , ρ ] k aa (cid:126) ¯ ω = ∂ α f k a (cid:126) ¯ ω (137) (cid:18) (cid:126) ¯ ω − ∆ (cid:15) k ◦ [ D α , ρ ] (cid:19) k ab = [ D α , ρ ] k ab (cid:126) ¯ ω − ∆ (cid:15) k ab = − i A α k ab ∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab (138)with a (cid:54) = b in the last equation.Expanding the first commutator in Eq. 135, (cid:20) D α , (cid:126) ¯ ω − ∆ (cid:15) k ◦ [ D α , ρ ] (cid:21) k aa = ∂ α ∂ α f k a (cid:126) ¯ ω + (cid:88) c A α k ac A α k ca ∆ f k ca (cid:126) ¯ ω + ∆ (cid:15) k ac + (cid:88) c A α k ca A α k ac ∆ f k ca (cid:126) ¯ ω − ∆ (cid:15) k ac (139)while the second commutator gives, for a (cid:54) = b , (cid:20) D α , (cid:126) ¯ ω − ∆ (cid:15) k ◦ [ D α , ρ ] (cid:21) k ab = − i (cid:126) ¯ ω A α k ab ( ∂ α ∆ f k ba ) − i A α k ab ; α ∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab − i A α k ab ( ∂ α ∆ f k ba ) (cid:126) ¯ ω − ∆ (cid:15) k ab − i A α k ab ( ∂ α ∆ (cid:15) k ab ) ∆ f k ba ( (cid:126) ¯ ω − ∆ (cid:15) k ab ) (140)Replacing Eqs. 139 and 140 in Eq. 135, σ βα α (¯ ω , ¯ ω ) = e (cid:126) (cid:126) ¯ ω ( (cid:126) ¯ ω + (cid:126) ¯ ω ) (cid:90) d d k (2 π ) d (cid:88) a ∂ β (cid:15) k a ∂ α ∂ α f k a + e (cid:126) (cid:126) ¯ ω + (cid:126) ¯ ω (cid:90) d d k (2 π ) d (cid:88) a,b A α k ba A α k ab ( ∂ β ∆ (cid:15) k ab )∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab − e (cid:126) (cid:126) ¯ ω (cid:90) d d k (2 π ) d (cid:88) a,b A β k ba A α k ab ∆ (cid:15) k ab ( ∂ α ∆ f k ba ) (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab − e (cid:126) (cid:90) d d k (2 π ) d (cid:88) a,b A β k ba A α k ab ∆ (cid:15) k ab ( ∂ α ∆ f k ba )( (cid:126) ¯ ω − ∆ (cid:15) k ab )( (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab ) − e (cid:126) (cid:90) d d k (2 π ) d (cid:88) a,b A β k ba A α k ab ; α ∆ f k ba ( (cid:126) ¯ ω − ∆ (cid:15) k ab )( (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab ) − e (cid:126) (cid:90) d d k (2 π ) d (cid:88) a,b A β k ba A α k ab ( ∂ α ∆ (cid:15) ab ) ∆ (cid:15) k ab ∆ f k ba ( (cid:126) ¯ ω − ∆ (cid:15) k ab ) ( (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab ) (141)Imposing time-reversal symmetry, we have (cid:15) − k a = (cid:15) k a and A − k ab = A k ba . Since it is only a function ofenergy, we have also f − k a = f k a . As a result, the first term in Eq. 141 vanishes due to the odd number ofderivatives in k , while the other terms are shown to have purely imaginary numerators by combining thecontributions at k and − k . For instance, 23 d d k (2 π ) d (cid:88) a,b A α k ba A α k ab ( ∂ β ∆ (cid:15) k ab )∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab = 12 (cid:90) d d k (2 π ) d (cid:88) a,b A α k ba A α k ab ( ∂ β ∆ (cid:15) k ab )∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab + 12 (cid:90) d d k (2 π ) d (cid:88) a,b A α k ba A α k ab ( ∂ β ∆ (cid:15) k ab )∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab = 12 (cid:90) d d k (2 π ) d (cid:88) a,b A α k ba A α k ab ( ∂ β ∆ (cid:15) k ab )∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab + 12 (cid:90) d d k (2 π ) d (cid:88) a,b A α − k ba A α − k ab ( − ∂ β ∆ (cid:15) − k ab )∆ f − k ba (cid:126) ¯ ω − ∆ (cid:15) − k ab = 12 (cid:90) d d k (2 π ) d (cid:88) a,b A α k ba A α k ab ( ∂ β ∆ (cid:15) k ab )∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab − (cid:90) d d k (2 π ) d (cid:88) a,b A α k ab A α k ba ( ∂ β ∆ (cid:15) k ab )∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab = (cid:90) d d k (2 π ) d (cid:88) a,b i Im {A α k ba A α k ab } ( ∂ β ∆ (cid:15) k ab )∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab (142)with the use of Eq. 32. This procedure is to be repeated for the remaining terms in Eq. 141.The next step is a partial fraction decomposition, σ βα α (¯ ω , ¯ ω ) = ie (cid:126) (cid:126) ¯ ω + (cid:126) ¯ ω (cid:90) d d k (2 π ) d (cid:88) a,b Im {A α k ba A α k ab } ( ∂ β ∆ (cid:15) k ab )∆ f k ba (cid:126) ¯ ω − ∆ (cid:15) k ab − ie (cid:126) (cid:126) ¯ ω (cid:90) d d k (2 π ) d (cid:88) a,b Im {A β k ba A α k ab } ( ∂ α ∆ f k ba ) (cid:18) (cid:126) ¯ ω + (cid:126) ¯ ω (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab − (cid:19) − ie (cid:126) (cid:90) d d k (2 π ) d (cid:88) a,b Im {A β k ba A α k ab } ( ∂ α ∆ f k ba ) (cid:18) (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω − ∆ (cid:15) k ab − (cid:126) ¯ ω + (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab (cid:19) − ie (cid:126) (cid:90) d d k (2 π ) d (cid:88) a,b Im {A β k ba A α k ab ; α } ∆ f k ba (cid:18) (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω − ∆ (cid:15) k ab − (cid:126) ¯ ω + (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab (cid:19) − ie (cid:126) (cid:90) d d k (2 π ) d (cid:88) a,b Im {A β k ba A α k ab } ( ∂ α ∆ (cid:15) ab ) ∆ f k ba × (cid:18) (cid:126) ¯ ω (cid:126) ¯ ω (cid:126) ¯ ω − ∆ (cid:15) k ab ) − (cid:126) ¯ ω + (cid:126) ¯ ω ( (cid:126) ¯ ω ) (cid:126) ¯ ω − ∆ (cid:15) k ab + (cid:126) ¯ ω + (cid:126) ¯ ω ( (cid:126) ¯ ω ) (cid:126) ¯ ω + (cid:126) ¯ ω − ∆ (cid:15) k ab (cid:19) (143)The second term in the second line of Eq. 143 can be singled out and identified as a Fermi surfacecontribution, ie (cid:126) (cid:126) ¯ ω (cid:90) d d k (2 π ) d (cid:88) a,b Im {A β k ba A α k ab } ( ∂ α ∆ f k ba ) = ie (cid:126) (cid:126) ¯ ω (cid:90) d d k (2 π ) d (cid:88) a,b i (cid:16) A β k ba A α k ab − A β k ab A α k ba (cid:17) ( ∂ α ∆ f k ba )= ie (cid:126) (cid:126) ¯ ω (cid:90) d d k (2 π ) d (cid:88) b ( − i ) (cid:2) A β , A α (cid:3) k bb ( ∂ α f k b ) = ie (cid:126) (cid:126) ¯ ω (cid:90) d d k (2 π ) d (cid:88) b F α β k b ( ∂ α f k b ) = ie (cid:126) (cid:126) ¯ ω F α α βB (144)with the use of Eqs. 32 and 102.Identifying the remaining integrals with the definitions in Eqs. 38-39 and 107, − (cid:126) i e σ βα α (¯ ω , ¯ ω ) = + ¯ ω ¯ ω Π βα α B (¯ ω ) − ¯ ω + ¯ ω ¯ ω Π βα α B (¯ ω + ¯ ω ) + ¯ ω + ¯ ω ¯ ω Π βα α B (¯ ω + ¯ ω ) − ω + ¯ ω Π α α β (¯ ω ) − ¯ ω + ¯ ω ¯ ω Π βα α (¯ ω ) − ¯ ω ¯ ω dd ¯ ω Π βα α (¯ ω ) + ¯ ω ¯ ω Π βα α (¯ ω )+ ¯ ω + ¯ ω ¯ ω Π βα α (¯ ω + ¯ ω ) − ¯ ω + ¯ ω ¯ ω Π βα α (¯ ω + ¯ ω ) (145)24nd using Eq. 106 to replace the derivative, we obtain (cid:126) i e σ βα α (¯ ω , ¯ ω ) = 1¯ ω + ¯ ω Π α α β (¯ ω ) + ¯ ω + ¯ ω ¯ ω Π βα α (¯ ω ) − ¯ ω ¯ ω Π α βα (¯ ω ) − ¯ ω + ¯ ω ¯ ω Π βα α (¯ ω + ¯ ω ) + ¯ ω + ¯ ω ¯ ω Π βα α (¯ ω + ¯ ω ) (146)in agreement with Eqs. 33-39. Appendix D: Integral evaluation for monolayer graphene
The evaluation of the integrals in Eqs. 29-30, 37-39 and 45-51 for monolayer graphene near the Dirac pointis relatively simple. Eqs 85 and 86 for the band energies and Eq. 87 for the non-Abelian Berry connectionare substituted and their derivatives (or generalized derivatives) are computed. Taking the relaxation-freelimit and considering T = 0, the I integrals that define the real part can be easily computed by convertingthe integral over the FBZ to polar coordinates and using the Dirac delta function in Eq. 54.It is important to note that the integrals must respect the same symmetries, as far as the tensor indices βα . . . α n are concerned, as the conductivities. This considerably reduces the number of integrals to compute,since, for graphene, σ xy ( ω ) = σ yx ( ω ) = 0 (147) σ xx ( ω ) = σ yy ( ω ) (148)at linear order. The second order is absent and at third order many tensor components vanish, σ xxxy ( ω , ω , ω ) = σ xxyx ( ω , ω , ω ) = σ xyxx ( ω , ω , ω ) = σ yxxx ( ω , ω , ω )= σ yyyx ( ω , ω , ω ) = σ yyxy ( ω , ω , ω ) = σ yxyy ( ω , ω , ω ) = σ xyyy ( ω , ω , ω ) = 0 (149)while the others are interrelated, σ xxxx ( ω , ω , ω ) = σ yyyy ( ω , ω , ω ) σ xxyy ( ω , ω , ω ) = σ yyxx ( ω , ω , ω ) (150) σ xyxy ( ω , ω , ω ) = σ yxyx ( ω , ω , ω ) σ xyyx ( ω , ω , ω ) = σ yxxy ( ω , ω , ω ) (151) σ xxxx ( ω , ω , ω ) = σ xxyy ( ω , ω , ω ) + σ xyxy ( ω , ω , ω ) + σ xyyx ( ω , ω , ω ) (152)with only three independent components remaining: σ xxyy , σ xyxy and σ xyyx . This leaves us with one andthree tensor elements to evaluate for each integral at linear and third order, respectively.For linear order, I xx ( ω ) = 14 π (cid:126) ω Θ( (cid:126) | ω | − | µ | ) (153)whose Hilbert transform H , H xx ( ω ) = − π − (cid:90) I xx ( ω (cid:48) ) ω (cid:48) − ω dω (cid:48) = 14 π (cid:126) ω ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12) (154)There is also the Fermi surface contribution F A , F xxA = − | µ | π (155)The case of the third order conductivity is entirely analogous, but there are, of course, more tensorcomponents to assess now, 25 xxyy ( ω ) = I xyxy ( ω ) = − (cid:126) v F π ( (cid:126) ω ) Θ( (cid:126) | ω | − | µ | ) (156) I xyyx ( ω ) = (cid:126) v F π ( (cid:126) ω ) Θ( (cid:126) | ω | − | µ | ) (157) I xxyy ( ω ) = I xyxy ( ω ) = I xyyx ( ω ) = (cid:126) v F π ( (cid:126) ω ) Θ( (cid:126) | ω | − | µ | ) (158) I xxyy ( ω ) = 3 (cid:126) v F π (cid:126) ω Θ( (cid:126) | ω | − | µ | ) (159) I xyxy ( ω ) = I xyyx ( ω ) = − (cid:126) v F π (cid:126) ω Θ( (cid:126) | ω | − | µ | ) (160) I xxyy ( ω ) = I xyxy ( ω ) = I xyyx ( ω ) = 0 (161) I xxyy ( ω ) = I xyxy ( ω ) = I xyyx ( ω ) = (cid:126) v F π ( (cid:126) ω ) Θ( (cid:126) | ω | − | µ | ) (162) I xxyy ( ω ) = − (cid:126) v F π ( (cid:126) ω ) Θ( (cid:126) | ω | − | µ | ) (163) I xyxy ( ω ) = I xyyx ( ω ) = (cid:126) v F π ( (cid:126) ω ) Θ( (cid:126) | ω | − | µ | ) (164)and the respective Hilbert transforms, H xxyy ( ω ) = H xyxy ( ω ) = − (cid:126) v F π ( (cid:126) ω ) ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12) − (cid:126) v F π | µ | (cid:126) ω (165) H xyyx ( ω ) = (cid:126) v F π ( (cid:126) ω ) ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12) + (cid:126) v F π | µ | (cid:126) ω (166) H xxyy ( ω ) = H xyxy ( ω ) = H xyyx ( ω ) = (cid:126) v F π ( (cid:126) ω ) ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12) + (cid:126) v F π | µ | (cid:126) ω (167) H xxyy ( ω ) = 3 (cid:126) v F π (cid:126) ω ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12) (168) H xyxy ( ω ) = H xyyx ( ω ) = − (cid:126) v F π (cid:126) ω ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12) (169) H xxyy ( ω ) = H xyxy ( ω ) = H xyyx ( ω ) = 0 (170) H xxyy ( ω ) = H xyxy ( ω ) = H xyyx ( ω ) = (cid:126) v F π ( (cid:126) ω ) ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12) + (cid:126) v F π | µ | ( (cid:126) ω ) (171) H xxyy ( ω ) = − (cid:126) v F π ( (cid:126) ω ) ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12) − (cid:126) v F π | µ | ( (cid:126) ω ) (172) H xyxy ( ω ) = H xyyx ( ω ) = (cid:126) v F π ( (cid:126) ω ) ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − | µ | (cid:126) ω + 2 | µ | (cid:12)(cid:12)(cid:12)(cid:12) + (cid:126) v F π | µ | ( (cid:126) ω ) (173)For the passage from the I to the H integrals, it is useful to hold in mind the following result for n > − π − (cid:90) Θ( ω (cid:48) − ∆ eff )( ω (cid:48) ) n ( ω − ω (cid:48) ) dω (cid:48) = 1 πω n ln (cid:12)(cid:12)(cid:12)(cid:12) (cid:126) ω − ∆ eff (cid:126) ω + ∆ eff (cid:12)(cid:12)(cid:12)(cid:12) + 2 π n − (cid:88) k =1 k −
1) ∆ k − eff ( (cid:126) ω ) n +1 − k (174)Finally, there is the Fermi surface contribution, F xxyyA = F xyxyA = F xyyxA = (cid:126) v F π | µ | (175)26eplacement of the integrals computed here in the Eqs. 41-44 of Section 3 reproduces the expression forthe third order conductivity derived by Cheng et al. [22]. References [1] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, Physical Review Letters (1961).[2] W. Kaiser and C. G. B. Garrett, Physical Review Letters (1961).[3] P. N. Butcher and D. Cotter, Elements of nonlinear optics (Cambridge university press, 1991).[4] P. N. Butcher and T. P. Mclean, Proceedings of the Physical Society (1964).[5] S. J. Lalama and A. F. Garito, Physical Review A (1979).[6] C. C. Teng and A. F. Garito, Physical Review Letters (1983).[7] P. N. Butcher and T. P. McLean, Proceedings of the Physical Society (1963).[8] C. Y. Fong and Y. R. Shen, Physical Review B (1975).[9] D. J. Moss, J. E. Sipe, and H. M. Van Driel, Physical Review B (1987).[10] D. J. Moss, E. Ghahramani, J. E. Sipe, and H. M. Van Driel, Physical Review B (1990).[11] J. E. Sipe and E. Ghahramani, Physical Review B (1993).[12] G. Ventura, D. Passos, J. Lopes Dos Santos, J. Viana Parente Lopes, and N. Peres, Physical Review B (2017).[13] D. Passos, G. Ventura, J. Viana Parente Lopes, J. Lopes Dos Santos, and N. Peres, Physical Review B (2018).[14] D. E. Parker, T. Morimoto, J. Orenstein, and J. E. Moore, Physical Review B (2019).[15] S. M. Jo˜ao and J. M. Lopes, arXiv (2018).[16] G. Ventura, D. Passos, J. Viana Parente Lopes, and J. Lopes Dos Santos, Journal of Physics CondensedMatter (2020).[17] C. Aversa and J. E. Sipe, Physical Review B (1995).[18] J. L. P. Hughes and J. E. Sipe, Physical Review B (1996).[19] J. E. Sipe and A. I. Shkrebtii, Physical Review B (2000).[20] C. Aversa, J. E. Sipe, M. Sheik-Bahae, and E. W. Van Stryland, Physical Review B (1994).[21] J. L. P. Hughes, Y. Wang, and J. E. Sipe, Physical Review B (1997).[22] J. L. Cheng, N. Vermeulen, and J. E. Sipe, New Journal of Physics (2014).[23] J. L. Cheng, N. Vermeulen, and J. E. Sipe, Physical Review B (2015).[24] E. I. Blount (Academic Press, 1962).[25] L. Allen and J. H. Eberly, Optical resonance and two-level atoms (Courier corporation, 1987).[26] T. Holder, D. Kaplan, and B. Yan, Physical Review Research (2020).[27] M. V. Berry, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences (1984). 2728] J. L. Cheng, J. E. Sipe, S. W. Wu, and C. Guo, APL Photonics (2019).[29] B. M. Fregoso, R. A. Muniz, and J. E. Sipe, Physical Review Letters (2018).[30] G. Ventura, D. Passos, J. Viana Parente Lopes, and J. Lopes dos Santos, arXiv (2021).[31] D. C. Hutchings, M. Sheik-Bahae, D. J. Hagan, and E. W. Van Stryland, Optical and Quantum Elec-tronics (1992).[32] M. Sheik-Bahae, D. J. Hagan, and E. W. Van Stryland, Physical Review Letters (1990).[33] M. Sheik-Bahae, D. C. Hutchings, D. J. Hagan, and E. W. Van Stryland, IEEE Journal of QuantumElectronics (1991).[34] W. R. Hannes and T. Meier, Physical Review B (2019).[35] A. H. Castro Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, Reviews of ModernPhysics (2009).[36] M. I. Katsnelson, Materials Today (2007).[37] S. A. Mikhailov, Europhysics Letters (2007).[38] S. A. Mikhailov and K. Ziegler, Journal of Physics Condensed Matter (2008).[39] E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, and S. A. Mikhailov, Physical Review Letters (2010).[40] R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. Peres, andA. K. Geim, Science (2008).[41] J. L. Cheng, J. E. Sipe, and C. Guo, Physical Review B (2019).[42] T. Jiang, D. Huang, J. Cheng, X. Fan, Z. Zhang, Y. Shan, Y. Yi, Y. Dai, L. Shi, K. Liu, C. Zeng, J. Zi,J. E. Sipe, Y. R. Shen, W. T. Liu, and S. Wu, Nature Photonics (2018).[43] H. M. Van Driel, J. E. Sipe, A. Hach´e, and R. Atanasov, Physica Status Solidi (B) Basic Research204