Nonlinear optical responses of crystalline systems: Results from a velocity gauge analysis
D. J. Passos, G. B. Ventura, J. M. Viana Parente Lopes, J. M. B. Lopes dos Santos, N. M. R. Peres
NNonlinear optical responses of crystalline systems:Results from a velocity gauge analysis
D. J. Passos, ∗ G. B. Ventura, J. M. Viana Parente Lopes, and J. M. B. Lopes dos Santos
Centro de F´ısica das Universidades do Minho e Porto andDepartamento de F´ısica e Astronomia, Faculdade de Ciˆencias,Universidade do Porto, 4169-007 Porto, Portugal
N. M. R. Peres
Centro de F´ısica das Universidades do Minho e Porto andDepartamento de F´ısica, Universidade do Minho, P-4710-057, Braga, Portugal
In this work, the difficulties inherent to perturbative calculations in the velocity gauge are ad-dressed. In particular, it is shown how calculations of nonlinear optical responses in the independentparticle approximation can be done to any order and for any finite band model. The procedure andadvantages of the velocity gauge in such calculations are described. The addition of a phenomeno-logical relaxation parameter is also discussed. As an illustration, the nonlinear optical response ofmonolayer graphene is numerically calculated using the velocity gauge.
I. INTRODUCTION
A common approach to understanding nonlinear opti-cal effects in atomic and condensed matter systems comesfrom a perturbation theory where the electric current isexpanded in powers of an external applied electric field ,assumed to be sufficiently weak for the expansion to bephysically meaningful. In this framework, one attemptsto calculate nonlinear optical response functions, usuallysecond or third order susceptibilities, the different fre-quency components of which describe a variety of physi-cal phenomena like the Kerr effect, harmonic generation,electro-optic effect, etc. The theory was first developed for atomic systems and,in the early nineties , extended to crystalline systems,characterized by electronic bands and the correspond-ing Bloch functions. The simpler calculations follow twoessential assumptions which will be adopted through-out this work: the independent particle approximation,where explicit electron-electron interactions are disre-garded (except possibly in the determination of the elec-tron bands, as in Density Functional Theory), and theelectric dipole approximation, where the response func-tions are taken to be local in space, a consequence of thelong wavelength limit of the applied electric field.Even in this simple approach, however, difficulties werefound regarding different ways to describe the perturba-tion. We can write the complete Hamiltonian of the crys-tal under the influence of the external field in two ways,either in the length gauge, H ( t ) = H ( r , p ) − q r · E ( t ) , (1)or in the velocity gauge or minimal coupling Hamiltonian, H ( t ) = H ( r , p − q A ( t )) , (2)where H is the unperturbed crystal Hamiltonian, q = − e is the charge of the electron and A ( t ) the vector potential, chosen so that E ( t ) = − ∂ t A ( t ). These twodescriptions can easily be shown to be equivalent andrelated by a time-dependent unitary transformation .The first attempts at computing nonlinear optical re-sponses in crystals made use of the minimal couplingmethod , since it has the advantage of retaining thetranslation symmetry of the Hamiltonian, only couplingstates with the same Bloch vector k . However, it pre-sented some serious difficulties as it seemed that responsefunctions computed in the velocity gauge diverged inthe limit of low frequencies, even for the case of zerotemperature insulators, where such divergences shouldclearly be absent. A more detailed analysis of the lin-ear response showed that these were only apparent diver-gences that could be removed by careful manipulationsand sum rules . As presented, these procedures werecumbersome and not easily generalizable to higher orderresponse functions.This led to a preference for the length gauge in non-linear optical response calculations . This gaugepresents its own difficulties, the most notorious one be-ing that the matrix elements of the position operator inthe Bloch basis are ill-defined until the thermodynamiclimit is taken and, even then, they can be defined onlyas a distribution. Inspired by Blount’s work , Aversa etal. pointed out that the position operator will appearin the calculations only inside commutators whose ma-trix elements are well defined and successfully developedthe theory in the length gauge. This formalism has sincebeen applied to various systems .More recently, and still within the length gauge ap-proach, simple expressions for the nonlinear conductivityto arbitrary order in the electric field were derived bythe present authors, making use of the covariant deriva-tive notation . These expressions are useful for inspec-tion and reduced the calculation procedure to a fairlystraightforward expansion of commutators. Although theidea of a “generalized derivative” was already around inthe literature , it involved separating the intra and inter- a r X i v : . [ c ond - m a t . o t h e r] M a r band components of the position operator and responsefunctions. The emphasis on these distinctions, originallymotivated by the intent of applying the equations to thespecial case of clean, cold semiconductors and to makeanalogies with atomic and free electron systems, madethe calculations less transparent, in our view.In this same work , similar and equivalent expressionswere derived using the velocity gauge, for the Hamilto-nian H ( r , p ) = p m + V L ( r ) , (3)where the second term is the periodic lattice potential.It was shown that, in the form derived from the velocitygauge, the expressions for the nonlinear optical coeffi-cients lose their validity if only a finite number of bandsaround the Fermi level are taken into account. This diffi-culty was recognized early on , and, together with theapparent infrared divergences, led to the velocity gaugebeing less adopted. The reasons for the failure of theseexpressions under a truncation in band space were alsounderstood and recently subjected to a more quan-titative analysis . One of the arguments consisted innoting that the sum rules that connected the expressionsin the two gauges seem to rely on commutator identitiessuch as (cid:2) r α , v β (cid:3) = i (cid:126) m δ αβ (4)( r and v , position and velocity operators), which requirean infinite number of bands to hold. This led to a com-mon misconception that the velocity gauge could only beproperly implemented if an infinite number of bands istaken into account . In fact, sum rules of greatergenerality have been constructed which remain true evenunder truncation of the bands . Nonetheless, various au-thors have indicated that the velocity gauge, if employed,would lead to different, unphysical, predictions .The fundamental difference between the two gaugesthat seems not to have been properly appreciated con-cerns the form of the perturbation. In the velocity gauge, the form of the perturbation depends explicitly on H ,unlike in the case of the length gauge. Expressions forthe nonlinear coefficients derived from the Hamiltonianin Eq. 3, cannot be truncated to a finite set of bands. Atruncation of H implies a different form of the perturba-tion, and leads to different expressions for the nonlinearconductivities in the velocity gauge, but not in the lengthgauge.In order to build the theory in its most general form,we will make no assumptions on the form of H , otherthan it has the periodicity of some Bravais lattice, sothat Bloch’s theorem applies and there is a well definedFirst Brillouin Zone (FBZ). The derived forms for thenonlinear conductivities will be completely equivalent to the ones obtained from the length gauge and can be ap-plied to any finite band model; no commutator identitiesof the kind of Eq. 4 are assumed.In the next section, the equivalence of the two gaugesis revisited, starting from an Hamiltonian with a finiteset of bands, and some formal concepts behind the ve-locity gauge formulation are reviewed. In section III, theequation of motion for the reduced density matrix in thevelocity gauge is introduced in a form more general thanpreviously presented. Recursively solving the equationof motion leads to the nonlinear conductivities. The ad-vantages and subtleties of using this gauge are discussedin Section IV, where it is argued that the velocity gaugeshould prove more efficient for numerical computations.In Section V, the introduction of a phenomenological re-laxation parameter is addressed, and proves to be lesstrivial than expected. As an example, in Section VIthe formalism is applied to third harmonic generationin graphene, in a full Brillouin zone calculation, and acomparison is made with already existing results in theliterature. Section VII contains some closing remarks. II. DENSITY MATRIX FORMALISM
In very general terms, which do not exclude finite bandmodels, the single particle unperturbed Hamiltonian H can be written as H = (cid:88) s (cid:90) d d k (2 π ) d | ψ k s (cid:105) [ H ] k ss (cid:104) ψ k s | , (5)where | ψ k s (cid:105) are the Bloch band states and [ H ] k ss (cid:48) ≡ (cid:15) k s δ ss (cid:48) , with (cid:15) k s the band energies.To represent the scalar potential, − q E · r , in the Blochbasis, we require from the start the infinite volume limit.We define the normalization of the Bloch wave functionsas (cid:104) ψ k (cid:48) s (cid:48) | ψ k s (cid:105) = (2 π ) d δ ss (cid:48) δ ( k − k (cid:48) ) . (6)where d is the dimensionality of the system. UsingBlount’s results for the matrix elements of the positionoperator , (cid:104) ψ k s | r | ψ k (cid:48) s (cid:48) (cid:105) = (2 π ) d [ δ ss (cid:48) ( − i ) ∇ k (cid:48) δ ( k (cid:48) − k )+ δ ( k − k (cid:48) ) ξ k ss (cid:48) ] , (7)where ξ k ss (cid:48) is the Berry connection , one obtains forthe single particle perturbed Hamiltonian in the lengthgauge (see Appendix B) H E = (cid:90) d d k (2 π ) d (cid:88) s,s (cid:48) | ψ k s (cid:105) [ δ ss (cid:48) (cid:15) k s − iq E ( t ) · D k ss (cid:48) ] (cid:104) ψ k s (cid:48) | (8)where the covariant derivative, D k ss (cid:48) , is defined by D k ss (cid:48) ≡ δ ss (cid:48) ∇ k − i ξ k ss (cid:48) . (9)In our previous paper , we discussed the equivalenceof the length and velocity gauges, starting from a the-ory with an infinite number of bands (Eq. 3); we showedthat a suitable truncation of final expressions for the non-linear optical response functions to a finite set of bandsleads to a reasonable approximation only in the lengthgauge. Therefore, if we want to formulate correctly avelocity gauge calculation with a finite set of bands, weshould take Eqs. 8 and 9 as our starting point, and ob-tain the single particle velocity gauge Hamiltonian, H A ,from a time dependent unitary transformation of H E . Inthis fashion, we preserve the equivalence of both descrip-tions, even when the sum over band indexes is finite. Inappendix B we derive H A = (cid:88) s,s (cid:48) (cid:90) d d k (2 π ) d | ψ k s (cid:105) H A k ss (cid:48) (cid:104) ψ k s (cid:48) | (10)with H A k ss (cid:48) = (cid:15) k s δ ss (cid:48) + V k ss (cid:48) ( t ) , (11)the time dependent perturbation, V k ss (cid:48) ( t ), being ex-pressed as a power series in the vector potential, V k ss (cid:48) ( t ) = ∞ (cid:88) n =1 ( − q ) n A α ( t ) . . . A α n ( t ) n ! h α ...α n k ss (cid:48) , (12)where h α ...α n k ss (cid:48) ≡ (cid:126) − n [ D α n k , [ . . . , [ D α k , H ]] ... ] ss (cid:48) . (13)An implicit summation over repeated Cartesian compo-nents α i is henceforth assumed. The coefficients in theexpansion are written explicitly in Eq. 13 in terms ofcommutators involving covariant derivatives.If we take the coefficient of the first order term in theexpansion of Eq. 13 as an example, h α k ss (cid:48) = (cid:126) − [ D α k , H ] ss (cid:48) = 1 (cid:126) ∂(cid:15) k s ∂k α δ ss (cid:48) − i (cid:126) ξ α k ss (cid:48) ( (cid:15) k s (cid:48) − (cid:15) k s )(14)These are the unperturbed velocity matrix elements,since v (0) k ss (cid:48) = − i (cid:126) − [ r , H ] k ss (cid:48) = (cid:126) − [ D k , H ] ss (cid:48) = h k ss (cid:48) . (15)Had we started from the Hamiltonian of Eq. 3, theperturbation expansion of Eq. 12 would be reduced tothe linear term in A ( t ); the second order term would have been a k independent constant, irrelevant for thedynamics, and all higher order terms would have beenzero . But to proceed in more general terms and ensureequivalence between length and velocity gauges for finiteband models, we must, for now, retain all the terms inthe series.For the analysis of the steady state response of the elec-tric current density, J , it is useful to rewrite the pertur-bation of Eq. 12 in frequency space, where the connectioncan already be made with the electric field componentsby E ( ω ) = iω A ( ω ), V k ss (cid:48) ( ω ) = ∞ (cid:88) n =1 (cid:90) + ∞−∞ dω π ... dω n π × ( iq ) n E α ( ω ) . . . E α n ( ω n ) n ! ω . . . ω n h α ...α n k ss (cid:48) × (2 π ) δ ( ω − ω − · · · − ω n ) (16)Having carefully defined the perturbation in the veloc-ity gauge, we can turn to the nonlinear optical responsefunctions, defined in frequency space by a similar expan-sion, (cid:104) J α (cid:105) ( ω ) = (cid:90) dω π σ (1) αβ ( ω ) E β ( ω ) (2 π ) δ ( ω − ω )+ (cid:90) dω π dω π σ (2) αβγ ( ω , ω ) E β ( ω ) E γ ( ω ) × (2 π ) δ ( ω − ω − ω )+ . . . (17)The ensemble average of the electric current density, asfor any observable, is given in terms of the density matrix ρ , (cid:104) J (cid:105) ( t ) = Tr[ J ρ ( t )] = q (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) v k ss (cid:48) Tr (cid:104) c † k s c k s (cid:48) ρ ( t ) (cid:105) = q (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) v k ss (cid:48) ρ k s (cid:48) s ( t ) (18)where the reduced density matrix (RDM) is defined as athe expectation of a product of a creation and a destruc-tion operator in Bloch states: ρ k s (cid:48) s ( t ) ≡ Tr (cid:104) c † k s c k s (cid:48) ρ ( t ) (cid:105) = (cid:104) c † k s c k s (cid:48) (cid:105) . (19)The standard density matrix formalism computes thenonlinear conductivities by expanding the RDM in pow-ers of the electric field. In the velocity gauge, however,the electric current is an explicitly time and field depen-dent observable . The velocity matrix elements are definedby − i (cid:126) − (cid:2) r , H A (cid:3) = (cid:126) − (cid:2) D , H A (cid:3) and also have to be ex-panded in powers of the electric field: v β k ss (cid:48) ( ω ) = ∞ (cid:88) n =0 (cid:90) + ∞−∞ dω π ... dω n π × ( iq ) n E α ( ω ) . . . E α n ( ω n ) n ! ω . . . ω n h βα ...α n k ss (cid:48) × (2 π ) δ ( ω − ω − · · · − ω n ) (20)The expansion must therefore be done in the density ma-trix and the velocity matrix elements simultaneously. Inthe absence of an external field, the current is (cid:104) J (cid:105) (0) = q (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) v (0) k ss (cid:48) ρ (0) k s (cid:48) s . (21)The first order response is (cid:104) J (cid:105) (1) = q (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) (cid:16) v (1) k ss (cid:48) ρ (0) k s (cid:48) s + v (0) k ss (cid:48) ρ (1) k s (cid:48) s (cid:17) , (22)and, in general, (cid:104) J (cid:105) ( n ) = q (cid:90) d d k (2 π ) d n (cid:88) p =0 (cid:88) ss (cid:48) v ( n − p ) k ss (cid:48) ρ ( p ) k s (cid:48) s (23)The expansion of the velocity matrix elements in theexternal field is already defined in Eq. 20. The expan-sion of the reduced density matrix involves solving itsequation of motion recursively. III. RECURSIVE SOLUTIONS TO THEEQUATION OF MOTION
In the absence of scattering, the equation of motion forthe reduced density matrix is i (cid:126) ∂ t ρ k ss (cid:48) = (cid:2) H A , ρ (cid:3) k ss (cid:48) (24)We can isolate the perturbation on the right hand side,( i (cid:126) ∂ t − ∆ (cid:15) k ss (cid:48) ) ρ k ss (cid:48) = [ V, ρ ] k ss (cid:48) (25)where ∆ (cid:15) k ss (cid:48) ≡ (cid:15) k s − (cid:15) k s (cid:48) .To solve the equation of motion perturbatively, it willbe written in frequency space, order by order in the elec-tric field, using the expansion of Eq. 16. For every order n , the reduced density matrix will be expressed in termsof its lower order terms. To alleviate notation and makethe recursion relation clearer, we factorize the electricfields and other common factors by defining, ρ ( n ) k ss (cid:48) ( ω ) ≡ (cid:90) dω π . . . dω n π ( iq ) n E α ( ω ) . . . E α n ( ω n ) ω . . . ω n × (2 π ) δ ( ω − ω − · · · − ω n ) ρ α ...α n k ss (cid:48) ( ω , . . . , ω n )(26)The goal is now to express the recursion relation be-tween objects of the form ρ α ...α n k ss (cid:48) ( ω , . . . , ω n ). In theabsence of a perturbation, the reduced density matrix issimply the Fermi-Dirac distribution, ρ (0) k ss (cid:48) = f ( (cid:15) k s ) δ ss (cid:48) = δ ss (cid:48) (cid:16) (cid:15) k s − µk B T (cid:17) , (27)which, when replaced in Eq. 21, implies (cid:104) J (cid:105) (0) = 0, asexpected.The first order term is ρ α k ss (cid:48) ( ω ) = (cid:104) h α k , ρ (0) k (cid:105) ss (cid:48) (cid:126) ω − ∆ (cid:15) k ss (cid:48) , (28)and the second order, ρ αβ k ss (cid:48) ( ω , ω ) = 1 (cid:126) ω + (cid:126) ω − ∆ (cid:15) k ss (cid:48) × (cid:18)(cid:104) h α k , ρ β k ( ω ) (cid:105) ss (cid:48) + 12 (cid:104) h αβ k , ρ (0) k (cid:105) ss (cid:48) (cid:19) (29)The pattern is already becoming clear. As an addi-tional example, the third order term has the form, ρ αβγ k ss (cid:48) ( ω , ω , ω ) = 1 (cid:126) ω + (cid:126) ω + (cid:126) ω − ∆ (cid:15) k ss (cid:48) × (cid:18)(cid:104) h α k , ρ βγ k ( ω , ω ) (cid:105) ss (cid:48) + 12 (cid:104) h αβ k , ρ γ k ( ω ) (cid:105) ss (cid:48) + 13! (cid:104) h αβγ k , ρ (0) k (cid:105) ss (cid:48) (cid:19) (30)Finally, to general order n , the perturbative solutionto the equation of motion is recursively expressed as ρ α ...α n k ss (cid:48) ( ω , . . . , ω n ) = 1 (cid:126) ω + · · · + (cid:126) ω n − ∆ (cid:15) k ss (cid:48) n (cid:88) m =1 m ! (cid:2) h α ...α m k , ρ α m +1 ...α n k ( ω m +1 , . . . , ω n ) (cid:3) ss (cid:48) (31)This recursion relation can be unfolded into lengthy ex-pressions and its structure analyzed in more detail. How-ever, we shall see that the real value of these expressionslies in their numerical evaluation, for which a recursion relation is sufficient.Applying Eq. 20 and Eq. 26 to the Eq. 23, the generalform of the nonlinear optical response functions in thevelocity gauge is obtained, σ ( n ) βα ...α n ( ω , . . . , ω n ) = i n q n +1 ω . . . ω n (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) n (cid:88) p =0 h βα ...α p k ss (cid:48) p ! ρ α p +1 ...α n k s (cid:48) s ( ω p +1 , . . . , ω n ) (32)This form of the nonlinear optical response functions willstill have to undergo the usual symmetrization procedureto ensure intrinsic permutation symmetry . Albeit triv-ial, this last step is a bit cumbersome to write down andwill be left implicit.These expressions are equivalent tothe ones we derived in a previous work using the lengthgauge. Although far more complicated, they have theiradvantages, which we will discuss in the next section.The equivalence of the results of the two approaches,length and velocity gauges, is guaranteed by their beingrelated by a unitary transformation (see Appendix B). Itcan also be explicitly shown by using very general sumrules to map the expressions for the nonlinear conductiv-ities in the velocity gauge to those of the length gauge,order by order, as shown in an appendix of our previouswork . The proof for first order responses is nonethelesspresented in appendix A, as an example of these sumrules, which generalize those derived in earlier works ,no longer rely on commutator identities (Eq. 4), and arevalid for a model with a finite number of bands. IV. LENGTH VS VELOCITY GAUGE
As usual, there are advantages and disadvantages as-sociated with any particular choice of gauge. By con-sidering the (exactly equivalent) forms of the nonlinearconductivities derived in the two gauges, the strengthsand weaknesses of each can be analyzed.A first look at Eq. 32 will immediately bring out theusual concerns with infrared divergences in the velocitygauge form, due to all the inverse frequency factors. We It is a consequence of the definition in Eq. 17 that only the sym-metric part of the response functions contributes to the integraland is therefore physical. emphasize again, however, that this expression is equiva-lent to the one obtained from the length gauge and there-fore these divergences are only apparent. Manipulatingthe expressions in the velocity gauge and using a series ofsum rules, it can be shown that Eq. 32 is divergence free.This approach was the one originally pursued , but thisuse of sum rules became rather pointless after the lengthgauge formulation had been developed . If the sum rulesare employed in the velocity gauge to remove apparentdivergences, one will simply arrive at the same expressionas obtained more straightforwardly in the length gauge.This equivalence through sum rules does not demand aninfinite number of bands, but it does put a constraint onthe use of Eq. 32, namely that the integration must bedone over the full FBZ to cancel divergences. An effec-tive continuum Hamiltonian describing a portion of theFBZ—like the Dirac Hamiltonian for graphene—, willnot suffice, since these sum rules rely on the periodicityin k space of the quantities involved.Having clarified this point, it can still be noted thatthe velocity gauge form is considerably more elaborate;less useful not only for inspection, but in an actual ana-lytical calculation. As an example, the expression of thesecond harmonic generation with all components alongthe x axis, in the length gauge , is , σ (2) xxx ( ω, ω ) = − q (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) h x k ss (cid:48) (cid:126) ω − ∆ (cid:15) k s (cid:48) s × (cid:20) D x k , (cid:126) ω − ∆ (cid:15) ◦ (cid:104) D x k , ρ (0) k (cid:105)(cid:21) s (cid:48) s (33)while in the velocity gauge, The symbol ◦ stands for the Hadamard product: ( A ◦ B ) ss (cid:48) = A ss (cid:48) B ss (cid:48) σ (2) xxx ( ω, ω ) = − q ω (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) (cid:16) h x k ss (cid:48) ρ xx k s (cid:48) s ( ω, ω )+ h xx k ss (cid:48) ρ x k s (cid:48) s ( ω ) + 12 h xxx k ss (cid:48) ρ (0) k s (cid:48) s (cid:19) (34)where we still have to write the density matrix compo-nents, ρ x k ss (cid:48) ( ω ) = (cid:104) h x k , ρ (0) k (cid:105) k ss (cid:48) (cid:126) ω − ∆ (cid:15) k ss (cid:48) (35) ρ xx k ss (cid:48) ( ω, ω ) = 12 (cid:126) ω − ∆ (cid:15) k ss (cid:48) × (cid:18) [ h x k , ρ x k ( ω )] k ss (cid:48) + 12 (cid:104) h xx k , ρ (0) k (cid:105) ss (cid:48) (cid:19) (36)This simple example illustrates that there is no advan-tage in doing the analytical calculations in the velocitygauge, although inspection of previous equations showsan interesting point: there are only simple poles in thevelocity gauge form ( (cid:126) ω − ∆ (cid:15) ) − , while in the lengthgauge, by differentiation, higher order poles emerge.Still, for analytical calculations, we would advocate themore transparent and simpler length gauge approach .The strength of the velocity gauge form lies in the dif-ferent arrangement of the commutators. The covariantderivatives are no longer applied to the density matrixin its recursion relation . Instead, they operate only onthe unperturbed Hamiltonian H in the determinationof the functions h k ss (cid:48) (Eq. 13), which are independent offrequency, temperature and chemical potential.A careful look at the algorithm of the previous sec-tion, shows that for the nonlinear conductivity of order n , there are n + 1 such functions to compute by suc-cessively applying a covariant derivative: h α ...α m k ss (cid:48) with m = 1 , ..., n + 1. In the previous example of second har-monic generation, these would be h x k ss (cid:48) , h xx k ss (cid:48) and h xxx k ss (cid:48) .Further reducing this algorithm to its fundamental in-gredients, we note that these calculations demand only aknowledge of two objects, which fully define the systemunder consideration: the dispersion relation (cid:15) k s and theBerry connection ξ k ss (cid:48) .Once these h α ...α n k ss (cid:48) functions are analytically deter-mined, the integrand in Eq. 32 or in Eq. 34 can be nu-merically evaluated at each k value, independently andquite easily. In fact, the procedure involves evaluatingthe analytic h α ...α n k ss (cid:48) functions and the Fermi-Dirac dis-tribution at the k point and then computing simple com-mutators and traces of numeric matrices. It has no nu-merical derivatives at all. This is in contrast with the In this aspect, the approach here has similarities with the oneemployed in . length gauge, where either the full expression of the re-sponse function is analytically calculated or numericalderivatives have to be applied in each step of the densitymatrix recursion relation. Either way, via the productrule and higher order poles, the number of complicatedterms to evaluate grows very fast with n in the lengthgauge approach.For this reason, the form of the nonlinear conductivityin Eq. 32, derived from the velocity gauge, should providea more efficient numerical approach. The authors haveimplemented numerically the expressions in both gaugesand done calculations on the nonlinear conductivity ofmonolayer graphene and observed that the computationtimes were indeed considerably smaller when Eq. 32 wasused. The velocity gauge results are presented in Sec-tion VI. V. PHENOMENOLOGICAL RELAXATIONPARAMETERS
In Eq. 32, like in all previous equations, the addition ofthe usual infinitesimal imaginary part to the frequencies, ω → ω + i + , is always implicit, as imposed by causal-ity. This provide us with well defined relaxation freeexpressions. From the numerical standpoint, the imagi-nary part must always be finite, but can be taken to besmaller than any other energy scale in the problem.However, one is also interested in modeling relaxationprocesses due to scattering of electrons with impurities,phonons and other electrons. A simple phenomenologicalapproach is to add one or more relaxation parameters tothe frequency poles. A common justification for the addi-tion of a relaxation parameter γ comes from consideringa scattering term in equation of motion , i (cid:126) ∂ t ρ k ss (cid:48) = [ H, ρ ] k ss (cid:48) − iγ ( ρ k ss (cid:48) − ρ eq k ss (cid:48) ) (37)If the perturbation is turned off, the density matrixrelaxes to the equilibrium distribution ρ eq . In the lengthgauge approach, ρ eq = ρ (0) . A simple rearrangement ,( i (cid:126) ∂ t + iγ − ∆ (cid:15) k ss (cid:48) ) ρ k ss (cid:48) = [ V, ρ ] k ss (cid:48) + iγρ (0) k ss (cid:48) (38)makes clear how this will impact the nonlinear conduc-tivities. The second term on the right hand side of Eq. 38is not relevant to any order n (cid:54) = 0, implying that the onlyalteration will be to add an imaginary constant to eachpole in frequency space, i (cid:126) ∂ t + iγ − ∆ (cid:15) k ss (cid:48) → (cid:126) ω + iγ − ∆ (cid:15) k ss (cid:48) (39) V here is not the same as in Eq. 20. It stands for the perturbationin the length gauge: V k ss (cid:48) ( t ) = − iqE α ( t ) D α k ss (cid:48) Again, using the second harmonic generation as an ex-ample, the scattering term will induce the following sim-ple modification of the expression in Eq. 33, σ (2) xxx ( ω, ω ) = − q (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) h x k ss (cid:48) (cid:126) ω + iγ − ∆ (cid:15) k s (cid:48) s × (cid:20) D x k , (cid:126) ω + iγ − ∆ (cid:15) ◦ (cid:104) D x k , ρ (0) k (cid:105)(cid:21) s (cid:48) s (40)In the case of the velocity gauge, this change is muchmore complicated, since the equilibrium distribution de-pends on the perturbation. Having in mind the connec-tion between the two gauges, it is easy to see that theequivalent ρ eq in the velocity gauge equation of motionshould be obtained by an unitary transformation of theFermi-Dirac distribution . This again translates into anexpansion in the vector potential, ρ eq k ss (cid:48) ( t ) = ∞ (cid:88) n =0 ( − q ) n A α ( t ) . . . A α n ( t ) n ! (cid:126) n × (cid:104) D α n k , (cid:104) . . . , (cid:104) D α k , ρ (0) k (cid:105)(cid:105) ... (cid:105) ss (cid:48) (41)If this equilibrium distribution is used, the results ob-tained will, indeed, be equivalent to the ones from thelength gauge. Nevertheless, as discussed in the previoussection, the main advantage of using the velocity gaugecomes from the absence of derivatives acting on the den-sity matrix, and once the term in Eq. 41 is added tothe equation of motion, that advantage, and the greaternumerical efficiency of this approach, are lost.As a consequence, we will follow a different phe-nomenology, more appropriate for our computations. Toeach frequency ω j we shall add a constant imaginarypart : (cid:126) ω j → (cid:126) ω j + iγ . This method comes naturallyfrom considering the adiabatic switching of the externalfields. It looks similar to the previous method, but theexpression for second harmonic generation (as would beobtained in the length gauge) σ (2) xxx ( ω, ω ) = − q (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) h x k ss (cid:48) (cid:126) ω + 2 iγ − ∆ (cid:15) k s (cid:48) s × (cid:20) D x k , (cid:126) ω + iγ − ∆ (cid:15) k ◦ (cid:104) D x k , ρ (0) (cid:105)(cid:21) s (cid:48) s (42)has now a factor of two in the relaxation parameter thatappears in one of the poles. This might seem like a slight The use of the Fermi-Dirac distribution as ρ eq in the velocitygauge leads to erroneous results, as already noted in , such asthe appearance of a nonzero electric current in the absence ofan applied electric field, by means of choosing a constant vectorpotential. This corresponds to choosing Γ ( n ) = nγ in the relaxation pa-rameters in Cheng et. al. . difference, but is a distinct phenomenology and, as weshall see, even for low values of γ it can lead to completelydifferent results for very low frequencies (cid:126) ω (cid:28) γ and ina region of width γ around resonances. VI. NUMERICAL IMPLEMENTATION: THIRDHARMONIC GENERATION IN GRAPHENE
In this section, the velocity gauge approach is testednumerically, by computing the third order conductivity ofa material whose nonlinear optical properties have beenthe subject of intensive research in recent years: mono-layer graphene .Monolayer graphene is a two dimensional sheet of car-bon atoms, displayed in an honeycomb lattice. It has twoatoms per Bravais lattice site. Here, we shall consideronly the simplest nearest neighbor tight binding modelthat describes the electronic properties of graphene , H = t (cid:20) φ ( k ) φ ∗ ( k ) 0 (cid:21) (43)where φ ( k ) = 1 + e − i k · a + e − i k · a = | φ ( k ) | e − iθ ( k ) (44)with t being the hopping parameter and a =(1 / , √ / a and a = ( − / , √ / a the basis vectorsof the Bravais lattice. From this model Hamiltonian, thedispersion relation and the Berry connection can be com-puted, (cid:15) k s = s t | φ ( k ) | (45) ξ k ss (cid:48) = − ss (cid:48) ∇ k θ (46)with band index s = ± − .The knowledge of the dispersion relation and theBerry connection is sufficient for a calculation of thenonlinear conductivity, independently of the gauge. Itcan be shown that a single independent component σ (3) xxxx ( ω, ω, ω ) describes third harmonic generation. Weshall confine ourselves to the study of this frequency com-ponent of the third order nonlinear conductivity.Results obtained from the standard length gauge ap-proach are presented in Fig. 1, for frequencies near the Under the approximation of no overlap between the Wannier or-bitals. Also, an additional constant term in the Berry connectionis neglected, since it is not relevant for frequencies near the Diracpoint.
ReIm0.00 0.05 0.10 0.15 0.20- 10- 5051015 ℏω / t σ xxxx ( ω , ω , ω ) / σ o (a) ReIm0.00 0.05 0.10 0.15 0.20- 10- 5051015 ℏω / t σ xxxx ( ω , ω , ω ) / σ o (b) FIG. 1: Frequency dependence of the third order non-linear conductivity of graphene, normalized to σ =3 q a t / π (cid:126) µ (same normalization used in ), at zero tem-perature. The parameters used were: (a) µ/t = 0 . γ/t =0 . µ/t = 0 . γ/t = 0 . γ introduced via a scattering term in the equationof motion (Eq. 38). ReIm0.064 0.065 0.066 0.067 0.068 0.069 0.070 0.071 - ℏω / t σ xxxx ( ω , ω , ω ) / σ o FIG. 2: Anomalous feature of the third order nonlinear con-ductivity of graphene from Fig. 1(b), in a region near thethree photon resonance 3 (cid:126) ω = 2 µ . Dirac point. In this case, the analytical form of the thirdorder conductivity in the Dirac point approximation wascalculated and then plotted. Also, γ is introduced viathe additional scattering term in the equation of motion(first type of phenomenological treatment described inthe previous section), as in Mikhailov’s work .These results are in agreement with analogous calcula-tions already in the literature , with the strongestfeature present at the three photon resonance 3 (cid:126) ω = 2 µ .In particular, Fig. 1(a) is identical to Fig. 3(b) andFig. 5 in refs and , respectively. Fig. 1(b) shows avery anomalous behavior at the resonance, which is al-ways present in a region | (cid:126) ω − µ | ≤ γ for small butfinite γ . Fig. 2 shows a close-up of this region. Thisstrange feature near resonance is analyzed in detail in ,where it is regarded as a prominent feature of potentialinterest, since in practice γ is always finite. However,a more careful analysis shows that despite the unusualshape of this feature, in the limit γ → + , the scatteringfree result in is indeed obtained (see Appendix C).To do the same calculations with the velocity gaugeapproach developed here, we evaluate analytically notthe full third order conductivity but only the followingfunctions h x k ss (cid:48) = a sin (cid:0) k x a (cid:1) C ss (cid:48) (cid:126) (cid:114) k x a ) + 4 cos (cid:0) k x a (cid:1) cos (cid:16) √ k y a (cid:17) (47) h xx k ss (cid:48) = a cos (cid:0) k x a (cid:1) C ss (cid:48) (cid:126) (cid:114) k x a ) + 4 cos (cid:0) k x a (cid:1) cos (cid:16) √ k y a (cid:17) (48) h xxx k ss (cid:48) = − a (cid:126) h x k ss (cid:48) h xxxx k ss (cid:48) = − a (cid:126) h xx k ss (cid:48) (49) C ≡ t × (cid:0) k x a (cid:1) + cos (cid:16) √ k y a (cid:17) − i sin (cid:16) √ k y a (cid:17) i sin (cid:16) √ k y a (cid:17) − (cid:0) k x a (cid:1) + cos (cid:16) √ k y a (cid:17) (50)All these analytical functions follow from direct evalu-ation of the commutators (Eq. 13) using Eqs. 45 and 46.At this point, the expression in Eq. 32 for the nonlin-ear conductivity (with n = 3) can be evaluated by nu-merical integration over the full FBZ. The phenomenol-ogy adopted here is the one that follows from adiabaticswitching, for the reasons discussed in the previous sec-tion.The velocity gauge results are in Fig. 3. The curvesare markedly different than the ones obtained from the Re (Velocity gauge)Im(Velocity gauge)Re (Length gauge)Im(Length gauge)0.00 0.05 0.10 0.15 0.20- 10- 5051015 ℏω / t σ xxxx ( ω , ω , ω ) / σ o (a) Re (Velocity gauge)Im(Velocity gauge)Re (Length gauge)Im(Length gauge)0.00 0.05 0.10 0.15 0.20- 10- 5051015 ℏω / t σ xxxx ( ω , ω , ω ) / σ o (b) FIG. 3: Frequency dependence of the third order non-linear conductivity of graphene, normalized to σ =3 q a t / π (cid:126) µ (same normalization used in ), at zero tem-perature. The parameters used were: (a) µ/t = 0 . γ/t =0 . µ/t = 0 . γ/t = 0 . γ was introduced by adiabaticswitching, simply replacing (cid:126) ω → (cid:126) ω + iγ . length gauge, especially for large γ . This should camewith no surprise, since a different phenomenologicaltreatment is adopted. To prove that the difference be-tween the two curves is only due to the way the relaxationparameters are introduced, the length gauge calculationswere performed again, now with this second type phe-nomenological treatment (using the third order analogueof Eq. 42) and included in the plots of Fig. 3. The resultsobtained in the two gauges are identical.
Of course, as we take the relaxation free limit γ → + the results of both phenomenological treatments also co-incide. However, for finite γ it may lead to some consider-ably different behavior of the nonlinear conductivity. Inparticular, we note that in the phenomenology adoptedhere ( (cid:126) ω i → (cid:126) ω i + iγ ), those anomalous features seen inFig. 1 and discussed in are not present. Instead, a moreplausible smooth curve is seen at the resonances.As a final remark, we emphasize that the results pre-sented here for the velocity gauge involved a completetwo-band tight binding model and an integration over the full FBZ, not an effective Dirac Hamiltonian, as inthe case of the length gauge results. This has some con-sequences. First, it had been suggested that a possiblesource for the two order of magnitude discrepancy be-tween theoretical and experimental results in graphenecould be the use of an effective Hamiltonian . Theagreement displayed in Fig. 3 establishes that the Diracpoint approximation is valid for the range of frequenciesadopted in previous studies . Second, the use ofcomplete bands will allow us to proceed towards higherfrequencies and study spectral regions where the DiracHamiltonian does not give an accurate description. VII. CONCLUSIONS
In summary, a velocity gauge approach to calculationsof nonlinear optical conductivities was developed in thiswork, within the density matrix formalism. It was shownthat, contrary to common belief, the results are the sameas those from the length gauge, as demanded by gauge in-variance. No difficulties come from applying this velocitygauge formalism to finite band models.The velocity gauge provides an efficient algorithm forcomputing nonlinear conductivities. The dispersion re-lation and the Berry connection of a crystal are the pre-requisites and define the crystalline system under study.A series of commutators (Eq. 13) can then be analyt-ically computed. From these, the conductivity can benumerically evaluated for any frequencies, temperature,chemical potential, relaxation parameter and to any or-der, without using any numerical derivatives.Previously, studies of nonlinear conductivities ofteninvolved writing down analytical expressions for thefull conductivity. For third order, this already be-comes rather cumbersome. Analytically, it is only re-ally tractable for simple effective continuum Hamiltoni-ans (such as the Dirac Hamiltonian). The velocity gaugeapproach developed here does not have such complica-tions, is easily generalizable to higher orders and is alwaysapplied to complete bands. This gauge provide us witha simple yet powerful method to compute the nonlinearoptical response functions of crystalline systems.
Appendix A: Equivalence of linear responses
As an illustration of the general sum rules implicit ina velocity gauge treatment, the expression for the linearconductivity in the velocity gauge σ (1) αβ ( ω ) = iq ω (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) h β k ss (cid:48) (cid:104) h α k , ρ (0) k (cid:105) s (cid:48) s (cid:126) ω − ∆ (cid:15) k s (cid:48) s + h βα k ss (cid:48) ρ (0) k s (cid:48) s (A1)will be shown to be equivalent to the one obtained froma length gauge approach,0 σ (1) αβ ( ω ) = − iq (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) h β k ss (cid:48) (cid:104) D α k , ρ (0) k (cid:105) s (cid:48) s (cid:126) ω − ∆ (cid:15) k s (cid:48) s (A2)To begin, the Jacobi identity is used to move the covari-ant derivative to the density matrix (cid:126) (cid:104) h α k , ρ (0) k (cid:105) = (cid:104) [ D α k , H ] , ρ (0) k (cid:105) = (cid:104)(cid:104) D α k , ρ (0) k (cid:105) , H (cid:105) + (cid:104) D α k , (cid:104) H , ρ (0) k (cid:105)(cid:105) = (cid:104)(cid:104) D α k , ρ (0) k (cid:105) , H (cid:105) (A3)where in the last step we took into account that the com-mutator of two diagonal matrices is zero (cid:104) H , ρ (0) k (cid:105) = 0.This leads to (cid:104) h α k , ρ (0) k (cid:105) ss (cid:48) = (cid:126) − (cid:104) D α k , ρ (0) k (cid:105) ss (cid:48) ∆ (cid:15) k s (cid:48) s (A4)The first term in in parenthesis of Eq. A1, dropping the k index for simplicity, becomes h βss (cid:48) (cid:2) h α , ρ (0) (cid:3) s (cid:48) s (cid:126) ω − ∆ (cid:15) s (cid:48) s = h βss (cid:48) (cid:126) − (cid:2) D α , ρ (0) (cid:3) s (cid:48) s ∆ (cid:15) ss (cid:48) (cid:126) ω − ∆ (cid:15) s (cid:48) s = h βss (cid:48) (cid:126) − (cid:104) D α , ρ (0) (cid:105) s (cid:48) s − ω h βss (cid:48) (cid:2) D α , ρ (0) (cid:3) s (cid:48) s (cid:126) ω − ∆ (cid:15) s (cid:48) s (A5)The second term when replaced in Eq. A1 will give thelength gauge result in Eq. A2. The remaining contribu-tions must therefore be zero and form our sum rule, iq ω (cid:90) d d k (2 π ) d (cid:88) ss (cid:48) (cid:16) h β k ss (cid:48) (cid:126) − (cid:104) D α k , ρ (0) k (cid:105) s (cid:48) s + h βα k ss (cid:48) ρ (0) k s (cid:48) s (cid:17) = 0(A6)This can be further simplified through h βαss (cid:48) ≡ (cid:126) − (cid:2) D α , h β (cid:3) ss (cid:48) (A7)to (cid:88) ss (cid:48) (cid:16) h βss (cid:48) (cid:126) − (cid:104) D α , ρ (0) (cid:105) s (cid:48) s + h βαss (cid:48) ρ (0) s (cid:48) s (cid:17) (A8)= (cid:126) − (cid:88) s (cid:104) D α , h β ρ (0) (cid:105) ss (A9)leading to the form iq (cid:126) ω (cid:90) d d k (2 π ) d (cid:88) s (cid:104) D α k , h β k ρ (0) k (cid:105) ss = 0 (A10) which can be recognized as a particular case of the sumrules identified in the appendix A of . The commutatorwith D α can be broken into two pieces, one involving theBerry connection which is trivially zero (the trace of aproper commutator is always zero) and another involvinga conventional derivative, iq (cid:126) ω (cid:90) d d k (2 π ) d (cid:88) s (cid:48) ∂∂k α (cid:16) h β k ss (cid:48) ρ (0) k s (cid:48) s (cid:17) = 0 (A11)This condition is always true since the functions h and ρ (0) are periodic in reciprocal space. The sum rule (andtherefore the equivalence between the results in the twogauges) is therefore trivially satisfied as long as the inte-gral is performed over the full FBZ. Appendix B: Expansion of H A on the vectorpotential In our previous paper , we discussed in detail theequivalence of the length and velocity gauges in the con-text of the unperturbed Hamiltonian of Eq. 3. In this ap-pendix, we review this equivalence, using from the starta representation of the crystal Hamiltonian in terms of aset of bands that can be finite.The representation of the position operator in theBloch basis requires the thermodynamic limit. Wechoose the following normalization for the Bloch states (cid:104) ψ k (cid:48) s (cid:48) | ψ k s (cid:105) = (2 π ) d δ ( k − k (cid:48) ) δ ss (cid:48) (B1)with the corresponding resolution of the identity (cid:88) s (cid:90) d d k (2 π ) d | ψ k s (cid:105) (cid:104) ψ k s | = ˆ1 (B2)where the sum over s may be over a finite set of bands.The unperturbed crystal Hamiltonian is H = (cid:88) s (cid:48) ,s (cid:90) d d k (2 π ) d | ψ k s (cid:105) [ H ] k ss (cid:48) (cid:104) ψ k s (cid:48) | (B3)with [ H ] k ss (cid:48) = (cid:15) k s δ ss (cid:48) , and (cid:15) k s the band energies. Theperturbation in the length gauge involves the positionoperator, H E = H − q E ( t ) · r (B4)Using Blount’s result for r in the thermodynamic limit ,we showed in our previous paper that (cid:104) ψ k s | r | ψ (cid:105) = (cid:88) s (cid:48) i D k ss (cid:48) (cid:104) ψ k s (cid:48) | ψ (cid:105) (B5)where the covariant derivative is defined by Eq. 9.From Eq. B5, we can give the position operator thefollowing representation: r = (cid:88) s,s (cid:48) (cid:90) d d k (2 π ) d | ψ k s (cid:105) i D k ss (cid:48) (cid:104) ψ k s (cid:48) | . (B6)1On first inspection, this equation might suggest that thisoperator is diagonal in Bloch momentum; it is not be-cause of the presence of the gradient with respect to k .Any observable that can be written as a differential op-erator acting on the wave function in a continuous basiswill have a similar representation .The full single particle Hamiltonian in the length gaugeis, therefore, H E = (cid:88) s,s (cid:48) (cid:90) d d k (2 π ) d | ψ k s (cid:105) [ δ ss (cid:48) (cid:15) k s − iq E ( t ) · D k ss (cid:48) ] (cid:104) ψ k s (cid:48) | (B7)The velocity gauge is obtained by a time dependent uni-tary transformation, | ψ A ( t ) (cid:105) = S ( t ) | ψ E ( t ) (cid:105) (B8) S ( t ) = (cid:88) s,s (cid:48) (cid:90) d d k (2 π ) d | ψ k s (cid:48) (cid:105) (cid:104) e − q (cid:126) A ( t ) · D k (cid:105) s (cid:48) s (cid:104) ψ k s | . (B9)The time evolution in this gauge is (cid:104) ψ k s | i (cid:126) ddt | ψ A ( t ) (cid:105) = (cid:104) ψ k s | S ( t ) ˆ H E S † ( t ) | ψ A ( t ) (cid:105) + (cid:104) ψ k s | i (cid:126) dS ( t ) dt S † ( t ) | ψ A ( t ) (cid:105) . (B10)The second term is seen to cancel the scalar potentialterm in the Hamiltonian, (cid:104) ψ k s | i (cid:126) dS ( t ) dt S † ( t ) | ψ A ( t ) (cid:105) = − i (cid:88) s (cid:48) q d A dt · D k ss (cid:48) (cid:104) ψ k s (cid:48) | ψ A ( t ) (cid:105) = i (cid:88) s (cid:48) q E ( t ) · D k ss (cid:48) (cid:104) ψ k s (cid:48) | ψ A ( t ) (cid:105) (B11)The velocity gauge Hamiltonian therefore becomes H A = (cid:88) s,s (cid:48) (cid:90) d d k (2 π ) d | ψ k s (cid:105) H A k ss (cid:48) (cid:104) ψ k s (cid:48) | (B12)with H A k ss (cid:48) ≡ (cid:88) r,r (cid:48) (cid:104) e − q (cid:126) A ( t ) · D k (cid:105) sr [ H ] k rr (cid:48) (cid:104) e q (cid:126) A ( t ) · D k (cid:105) r (cid:48) s (cid:48) At this point, we make use of the following identity forany two operators ˆ B and ˆ C , e ˆ C ˆ B e − ˆ C = ˆ B + (cid:104) ˆ C, ˆ B (cid:105) + 12! (cid:104) ˆ C, (cid:104) ˆ C, ˆ B (cid:105)(cid:105) + . . . (B13)and apply it to Eq. B12 with ˆ B = H and ˆ C = − q (cid:126) A ( t ) · D k , to get H A k ss (cid:48) = ∞ (cid:88) n =0 ( − q ) n A α ( t ) . . . A α n ( t ) (cid:126) n n ! × [ D α n k , [ . . . , [ D α k , H ]] ... ] ss (cid:48) = (cid:15) k s δ ss (cid:48) + ∞ (cid:88) n =1 ( − q ) n A α ( t ) . . . A α n ( t ) (cid:126) n n ! × [ D α n k , [ . . . , [ D α k , H ]] ... ] ss (cid:48) (B14) Appendix C: Brief note on a phenomenologicalfeature In , Mikhailov pointed out that the feature observedin Fig. 2 was due mainly to a term of the form ∼ Im γ (3 (cid:126) ω − µ + iγ ) (C1)In the limit γ →
0, this should be considered a distribu-tion, to be integrated over frequency. One knows thatlim γ → Im 13 (cid:126) ω − µ + iγ (C2)corresponds to a Dirac delta function, so one may askto what distribution corresponds Eq. C1. Since we canwritelim γ → Im γ (3 (cid:126) ω − µ + iγ ) = − (cid:126) lim γ → Im ∂ ω γ (cid:126) ω − µ + iγ (C3)upon integration with a test function f ( ω ), we obtain acontribution ∝ lim γ → γ f (cid:48) (2 µ/ (cid:126) ) which is always zero. Inthis limit, the term in Eq. C1 has no weight at all. It willnot show up in an integration over ω . Similar remarks canbe made concerning the other resonances in the nonlinearconductivity, computed via the standard phenomenolog-ical approach; this is how the two phenomenologies dis-cussed in Section V yield different results for any finite γ, yet, nevertheless, become identical in the γ → + limit.2 The more familiar case is that of the momentum p x = − i (cid:126) (cid:82) dx | x (cid:105) ∂ x (cid:104) x | . ∗ Electronic address: [email protected] Yuen-Ron Shen. The principles of nonlinear optics.
NewYork, Wiley-Interscience, 1984, 575 p. , 1984. JE Sipe and Ed Ghahramani. Nonlinear optical responseof semiconductors in the independent-particle approxima-tion.
Physical Review B , 48(16):11705, 1993. Claudio Aversa and JE Sipe. Nonlinear optical suscep-tibilities of semiconductors: Results with a length-gaugeanalysis.
Physical Review B , 52(20):14636, 1995. GB Ventura, DJ Passos, JMB Lopes dos Santos,JM Viana Parente Lopes, and NMR Peres. Gauge covari-ances and nonlinear optical responses.
Physical Review B ,96(3):035431, 2017. VN Genkin and PM Mednis. Contribution to the theory ofnonlinear effects in crystals with account taken of partiallyfilled bands.
Sov. Phys. JETP , 27:609, 1968. James LP Hughes and JE Sipe. Calculation of second-order optical response in semiconductors.
Physical ReviewB , 53(16):10751, 1996. JE Sipe and AI Shkrebtii. Second-order optical responsein semiconductors.
Physical Review B , 61(8):5337, 2000. Ibraheem Al-Naib, JE Sipe, and Marc M Dignam. Highharmonic generation in undoped graphene: Interplayof inter-and intraband dynamics.
Physical Review B ,90(24):245423, 2014. F Hipolito, Thomas G Pedersen, and Vitor M Pereira.Nonlinear photocurrents in two-dimensional systems basedon graphene and boron nitride.
Physical Review B ,94(4):045434, 2016. EI Blount. Formalisms of band theory.
Solid state physics ,13:305–373, 1962. JL Cheng, N Vermeulen, and JE Sipe. Dc current inducedsecond order optical nonlinearity in graphene.
Optics ex-press , 22(13):15868–15876, 2014. Jin Luo Cheng, Nathalie Vermeulen, and JE Sipe. Third-order nonlinearity of graphene: Effects of phenomenologi-cal relaxation and finite temperature.
Physical Review B ,91(23):235320, 2015. Fred Nastos, Bernd Olejnik, Karlheinz Schwarz, andJE Sipe. Scissors implementation within length-gaugeformulations of the frequency-dependent nonlinear op-tical response of semiconductors.
Physical Review B ,72(4):045223, 2005. Claudio Aversa, JE Sipe, M Sheik-Bahae, andEW Van Stryland. Third-order optical nonlinearitiesin semiconductors: The two-band model.
Physical Review B , 50(24):18073, 1994. K Rzazewski and Robert W Boyd. Equivalence of inter-action hamiltonians in the electric dipole approximation.
Journal of modern optics , 51(8):1137–1147, 2004. Alireza Taghizadeh, F Hipolito, and TG Pedersen. Linearand nonlinear optical response of crystals using length andvelocity gauges: Effect of basis truncation. arXiv preprintarXiv:1710.01300 , 2017. Kuljit S Virk and JE Sipe. Semiconductor optics in lengthgauge: A general numerical approach.
Physical Review B ,76(3):035213, 2007. JL Cheng, Nathalie Vermeulen, and JE Sipe. Third orderoptical nonlinearity of graphene.
New Journal of Physics ,16(5):053014, 2014. Michael V Berry. Quantal phase factors accompanying adi-abatic changes. In
Proceedings of the Royal Society of Lon-don A: Mathematical, Physical and Engineering Sciences ,volume 392, pages 45–57. The Royal Society, 1984. Di Xiao, Ming-Che Chang, and Qian Niu. Berry phaseeffects on electronic properties.
Reviews of modern physics ,82(3):1959, 2010. Sergey A Mikhailov. Quantum theory of the third-ordernonlinear electrodynamic effects of graphene.
Physical Re-view B , 93(8):085403, 2016. JL Cheng, N Vermeulen, and JE Sipe. Second order opticalnonlinearity of graphene due to electric quadrupole andmagnetic dipole effects.
Scientific Reports , 7, 2017. MM Glazov and SD Ganichev. High frequency electric fieldinduced nonlinear effects in graphene.
Physics Reports ,535(3):101–138, 2014. SA Mikhailov. Non-linear electromagnetic response ofgraphene.
EPL (Europhysics Letters) , 79(2):27002, 2007. SA Mikhailov. Nonperturbative quasiclassical theory of thenonlinear electrodynamic response of graphene.
PhysicalReview B , 95(8):085432, 2017. A Marini, JD Cox, and FJ Garc´ıa de Abajo. Theoryof graphene saturable absorption.
Physical Review B ,95(12):125408, 2017. NMR Peres, Yu V Bludov, Jaime E Santos, Antti-PekkaJauho, and MI Vasilevskiy. Optical bistability of graphenein the terahertz range.
Physical Review B , 90(12):125425,2014. AH Castro Neto, F Guinea, Nuno MR Peres, Kostya SNovoselov, and Andre K Geim. The electronic propertiesof graphene.