Gauge covariances and nonlinear optical responses
G. B. Ventura, D. J. Passos, J. M. B. Lopes dos Santos, J. M. Viana Parente Lopes, N. M .R. Peres
aa r X i v : . [ c ond - m a t . o t h e r] J un Gauge covariances and nonlinear optical responses
G. B. Ventura, ∗ D. J. Passos, 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
J. M. Viana Parente Lopes
Centro de F´ısica das Universidades do Minho e PortoDepartamento de Engenharia F´ısica, Faculdade de Engenharia 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
The formalism of the reduced density matrix is pursued in both length and velocity gauges ofthe perturbation to the crystal Hamiltonian. The covariant derivative is introduced as a convenientrepresentation of the position operator. This allow us to write compact expressions for the reduceddensity matrix in any order of the perturbation which simplifies the calculations of nonlinear opti-cal responses; as an example, we compute the first and third order contributions of the monolayergraphene. Expressions obtained in both gauges share the same formal structure, allowing a com-parison of the effects of truncation to a finite set of bands. This truncation breaks the equivalencebetween the two approaches: its proper implementation can be done directly in the expressionsderived in the length gauge, but require a revision of the equations of motion of the reduced densitymatrix in the velocity gauge.
I. INTRODUCTION
The calculation of nonlinear optical (NLO) coeffi-cients in crystals has seen a renewed impetus, spurredby the strong nonlinear properties of layered materialslike graphene [1–8].Perturbative calculations of NLO coefficients in bulksemiconductors, with a full quantum treatment of mat-ter, date back to early nineties of the previous century,and have not been entirely trouble free [9, 10]. In thelong wavelength limit—in which the spatial dependenceof the radiation electric field is neglected—, there aretwo representations of the radiation field: by a timedependent vector potential A ( t ), with the electric fieldgiven by E ( t ) = − ∂ A /∂t ; by the electric dipole scalarpotential, V ( r ) = e E ( t ) · r . The advantage of the firstmethod, known as the velocity gauge, is that the per-turbation introduces no extra spatial dependence to thecrystal Hamiltonian, thus preserving the crystal’s trans-lational symmetry. This leads to a decoupling of thesystem’s response in momentum space: it becomes asum of independent contributions of each k value in theBrillouin zone. Early attempts to calculate NLO coef-ficients using this approach were, however, plagued byunphysical contributions, diverging at low frequencies[9]. Several authors addressed this issue by separatingthe treatment of inter and intra band contributions, us-ing time-dependent basis sets [10, 11]. Later Aversa and ∗ corresponding author: [email protected] Sipe [12] revisited the problem, emphasizing the gaugefreedom that allows you to choose either form of the cou-pling to the radiation field. They recognized that theunphysical divergences mentioned above actually havecoefficients that are exactly zero, expressing sum rulesthat they derived explicitly in first order response, andclaimed to hold in all orders. Because these sums rulesare easily violated in approximations, they end up advo-cating using the scalar potential method, also referredto as the length gauge, in actual calculations.Similar problems were found in earlier calculationsof NLO response of atoms [13]. As far back as 1951,ref. [14], W. Lamb Jr. recommended as more convenientthe scalar potential gauge in perturbative calculationsof the Hydrogen atom fine structure, and, for a while,the view that the two choices of gauge lead to differentresults was widely held [13].The obvious advantage of the scalar potential gauge isthat it is written in terms of a gauge invariant entity, theelectric field, even though it expresses a specific choiceof gauge for the electromagnetic field. But because thescalar potential contains the position variable, r , theperturbation is no longer diagonal in Bloch momentumspace, and couples different k values. Furthermore, theposition operator is highly singular in momentum space,and its matrix elements can only be properly defined inthe infinite crystal limit.This choice of representation was used on the recentreduced density matrix (RDM) calculations of the non-linear optical response of graphene [4, 6, 8]. In ref. [4],the derivative term in the RDM equations of motionwas removed by means of a k -space translation thus de-coupling them in crystal momentum space; this is notwithout cost, as the system’s response is now expressedin terms of both the A and the E fields. In their subse-quent work [6], the authors retained this derivative termas well as introduced relaxation terms to the RDM equa-tions of motion. A different approach was proposed byMikhailov [8], who avoided the problem of the singu-lar intra band term of the position operator by usinga finite wavelength perturbation that satisfies periodicboundary conditions, V ( r ) = − e (cid:0) V q e i q · r + V ∗ q e − i q · r (cid:1) , (1)and taking the limit q → ψ k s → e iθ s ( k ) ψ k s , since any expectation value is necessarily independentfrom the choice of θ s ( k ). This sets the transformationlaw of any observable’s matrix elements to A kk ′ ss ′ → e − i ( θ s ( k ) − θ s ′ ( k ′ )) A kk ′ ss ′ . Striving to make this property explicit leads to the con-cept of the covariant derivative, which will be used toderive a considerably simpler form of the RDM formal-ism in the length gauge.This paper is organized as follows. In Section II, wepresent an overview of some concepts on gauge invari-ance and the key ideas regarding the crystal Hamilto-nian. We then introduce the covariant derivative in mo-mentum space, which captures the phase freedom men-tioned in the previous paragraph. Section III is dedi-cated to the reduced density matrix. We shall derivethe RDM equations of motions and, more importantly,the relation between these two objects. The rest of thesection is dedicated to the equivalence between observ-ables in the two formalisms. In Section IV, we writethe solutions to the RDM equations of motion. As aproof of concept, Section V is dedicated to the studyof graphene’s current response by means of the scalarpotential formalism [4, 6, 8]. This is followed by Sec-tion VI, where we discuss the breakdown of the scalarpotential/vector potential equivalence, upon truncationof the expressions for the current for a finite set of bands[15]. The last section is dedicated to a summary of ourresults.
II. A GENERAL DESCRIPTION
The two possible representations of the uniform elec-tric field entail two different, but equivalent, ways towrite the many-body Hamiltonian. One can either addthe dipole interaction to the single particle Hamiltonian of the unperturbed system, H , [16] H E ( t ) = Z d d r Ψ † ( r ) (cid:20) H (cid:18) r , ∇ i (cid:19) + e E ( t ) · r (cid:21) Ψ( r ) , (2)or use the minimum coupling procedure in H , H A ( t ) = Z d d r Ψ † ( r ) (cid:20) H (cid:18) r , ∇ i + e ~ A ( t ) (cid:19)(cid:21) Ψ( r ) . (3)The electron field, Ψ( r ), and its Hermitian conjugate,Ψ † ( r ), satisfy the usual anti-commutation relations.The many-body state vector in the vector potentialapproach, | ψ ( t ) i , evolves in time according to the Hamil-tonian H A ( t ) i ~ ∂ | ψ ( t ) i ∂t = H A ( t ) | ψ ( t ) i . (4)A second state vector (cid:12)(cid:12) ¯ ψ ( t ) (cid:11) , obtained via a time-dependent unitary transformation of | ψ ( t ) i , (cid:12)(cid:12) ¯ ψ ( t ) (cid:11) = U ( t ) | ψ ( t ) i , (5)has an equation of motion i ~ ∂ (cid:12)(cid:12) ¯ ψ ( t ) (cid:11) ∂t = (cid:20) U ( t ) H A ( t ) U † ( t ) + i ~ d U ( t ) dt U † ( t ) (cid:21) (cid:12)(cid:12) ¯ ψ ( t ) (cid:11) . If the unitary transformation is chosen as U ( t ) = exp (cid:20) i e ~ Z d d r A ( t ) · r ρ ( r ) (cid:21) , (6)where ρ ( r ) := Ψ † ( r )Ψ( r ) is the density operator, it isstraightforward to show that U ( t ) H A ( t ) U † ( t ) + i ~ d U ( t ) dt U † ( t ) = H E ( t ) , (7)implying that (cid:12)(cid:12) ¯ ψ ( t ) (cid:11) is the state vector in scalar poten-tial gauge.Observables in the two gauges are also related by aunitary transformation, O E := U ( t ) O A ( t ) U † ( t ). Theexception is the Hamiltonian; because the unitary trans-formation, U ( t ), is time dependent, the time evolu-tion operator in the length gauge, H E ( t ), is not sim-ply U ( t ) H A ( t ) U † ( t ), but has an additional term involv-ing the time derivative of U ( t ), Eq. (7). The existenceof this transformation between the two descriptions ofthe radiation field establishes their complete equivalence[12, 13].Next we recall some important results of electroneigenstates in an unperturbed crystal. The single parti-cle Schr¨odinger equation is [17] H ψ k s ( r ) = ǫ k s ψ k s ( r ) , (8)with H = ~ m (cid:18) ∇ i (cid:19) + V ( r ) , (9)and V ( r ) = V ( r + R ), for R any Bravais lattice vector.According to Bloch’s theorem, the eigenfunctions havethe form of a plane wave times a periodic function, ψ k s ( r ) = e i k · r u k s ( r ) , (10)allowing the eigenvalue problem to be expressed in termsof the k -dependent Hamiltonian, H ( k ) := e − i k · r H e i k · r , H ( k ) u k s ( r ) = " ~ m (cid:18) ∇ i + k (cid:19) + V ( r ) u k s ( r )= ǫ k s u k s ( r ) . (11)The function u k s ( r ) is a periodic function in the realspace unit cell, u k s ( r ) = u k s ( r + R ) . (12)Each k -point in the First Brillouin Zone (FBZ) definesan Hamiltonian operator, H ( k ), that acts on functionswhose domain is the real space unit cell, and that satisfythe boundary condition of Eq. (12). The eigenfunctionsof H k , {| u k s i , s = 0 , , . . . } are a basis of such func-tions. We can assume this basis to be orthonormal withan inner product defined as the integral over the realspace unit cell (volume v c ), h u k s | u k s ′ i = 1 v c Z uc d d r u ∗ k s ( r ) u k s ′ ( r ) = δ ss ′ . (13)Different values of k have different basis, for they areeigenfunctions of different Hamiltonians. Here, theBloch wave vector is a continuous parameter, even in thefinite volume crystal, as the eigenvalues in Eq. (11) arewell defined for every k in the FBZ. The k value selectionby periodic boundary conditions only involves the planewave factor of the Bloch function, and has no bearing onthe periodic part, u k s ( r ) [17]. As such, derivatives withrespect to k of u k s ( r ) are always well defined whereasderivatives of plane wave factors require the infinite vol-ume limit. We shall work in this limit from the start,due to the difficulties of properly defining the positionoperator, r , in a finite system with periodic boundaryconditions.In the Ω → ∞ limit, (Ω, the volume of the crystal)the momentum sums are replaced by d -dimensional in-tegrals over the FBZ. The many-body crystalline Hamil-tonian then reads as H = X s Z d d k (2 π ) d ǫ k s c † k s c k s , (14)for c † k s , c k s the creation and destruction operators ofBloch states, c † k s = Z d d r ψ k s ( r ) Ψ † ( r ) . (15)The Bloch state orthogonality relation, the anti-commutation relations, and the lattice sum rule are suit-ably modified, h ψ k s | ψ k ′ s ′ i = (2 π ) d δ ss ′ δ ( k − k ′ ) , (16) (cid:8) c k s , c † k ′ s ′ (cid:9) = (2 π ) d δ ss ′ δ ( k − k ′ ) , (17) X R e i ( k − k ′ ) · R = (2 π ) v C δ ( k − k ′ ) , (18)so that the operator c † k s c k s is a density in momentumspace and not a dimensionless number operator as inthe finite volume case.In the scalar potential approach, the perturbation iswritten in terms of the position operator, r . Its matrixelements are ill-defined in the finite volume system, butcan be computed for Ω → ∞ . In that limit they read[18], r kk ′ ,ss ′ = δ ss ′ (2 π ) d ( − i ) ∇ k ′ δ ( k ′ − k )+(2 π ) d δ ( k ′ − k ) ξ k ′ ss ′ , (19)where the Berry connection, ξ k ss ′ , is defined as a scalarproduct in the real space unit cell, independent of thecrystal’s volume, ξ k ss ′ := i h u k s |∇ k u k s ′ i (20)= iv C Z uc d d r u ∗ k s ( r ) ∇ k u k s ′ ( r ) . (21)The somewhat awkward looking expression of Eq. (19),can be cast in a more transparent form if we bear inmind that, for a continuous non-normalizable basis, thematrix elements of an operator are a kernel of an integraltransform. In the Bloch representation, a general singleparticle state is represented asΨ( r ) = X s Z d d k (2 π ) d Φ s ( k ) ψ k s ( r ) , for Φ s ( k ) = h ψ k s | Ψ i . The wave function for the state r | Ψ i is h ψ k s | r | Ψ i = X s ′ Z d d k ′ (2 π ) d r kk ′ ,ss ′ Φ s ′ ( k ′ )= i X s ′ ( δ ss ′ ∇ k − i ξ k ss ′ ) Φ s ′ ( k ) , where, to reach the final expression, we have integratedby parts. The integration is over the FBZ, and, giventhe periodicity of any function of k , φ ( k ) = φ ( k + G ),there are no surface terms. This prompts us do definethe covariant derivative operator [6], D k ss ′ := δ ss ′ ∇ k − i ξ k ss ′ , (22)such that h ψ k s | r | Ψ i = X s ′ i D k ss ′ Φ s ′ ( k ) , (23) i.e., the position operator is i D k ss ′ in the Bloch repre-sentation [19]. The designation of covariant refers to itsbehavior under a local gauge transformation in momen-tum space, u k s → e iθ s ( k ) u k s , (24)for which D k ss ′ → ˜ D k ss ′ := e − iθ s ( k ) D k ss ′ e iθ s ′ ( k ) . = e − i ( θ s ( k ) − θ s ′ ( k )) D k ss ′ . The gradient term of the phase θ s ′ ( k ) is canceled bythe transformation of the Berry connection; this is incomplete parallel to the definition of covariant derivativein gauge theories in real space.In the vector potential approach, the perturbationis written in terms of the velocity matrix elements, v kk ′ ss ′ := p kk ′ ss ′ /m e , which are diagonal in k -space,and expressible as matrix elements in the basis of peri-odic functions, v kk ′ ss ′ = (2 π ) d δ ( k − k ′ ) v k ss ′ , (25) v k ss ′ = ~ m e h u k s | ( − i ∇ + k ) | u k s ′ i , (26)Since the velocity operator is quite generally( i ~ ) − [ r , H ], it is not surprising to find that v k ss ′ = 1 ~ [ D k , H ( k )] ss ′ , (27)= 1 ~ [ δ ss ′ ∇ k ǫ k s − i ( ǫ k s ′ − ǫ k s ) ξ k ss ′ ] , (28)This turns out to be the expression of a more generalresult, which will prove useful later and which we nowdiscuss.Any matrix in the space generated by the basis {| u k s i , s = 1 , , . . . } , parametrized by k , defines an op-erator in band space by O ( k ) = X ss ′ | u k s i O k ss ′ h u k s ′ | , (29) O k ss ′ = h u k s | O ( k ) | u k s ′ i . (30)Two examples of this are the Hamiltonian, H ( k ) := ǫ k s δ ss ′ and the velocity operator. Consider the matrixelements of the k -derivative of O ( k ),[ ∇ k O ( k )] ss ′ := h u k s | ∇ k O ( k ) | u k s ′ i , = ∇ k h u k s | O ( k ) | u k s ′ i− h∇ k u k s | O ( k ) | u k s ′ i− h u k s | O ( k ) |∇ k u k s ′ i . (31) The first term on the right hand side is simply the k -gradient of the matrix element, ∇ k O k ss ′ . The othertwo can be expressed in terms of the Berry con-nection, by application of the completeness relations P r | u k r i h u k r | = ˆ1, h∇ k u k s | O ( k ) | u k s ′ i = i X r ξ k sr O k rs ′ , (32) h u k s | O ( k ) |∇ k u k s ′ i = − i X r O k sr ξ k rs ′ . (33)The matrix element, Eq. (31), then reads as a commu-tator with the covariant derivative,[ ∇ k O ( k )] ss ′ = ∇ k O k ss ′ − i [ ξ k , O ( k )] ss ′ , (34)= (cid:2) D k , O ( k ) (cid:3) ss ′ . (35)which can be alternatively represented in operator form, ∇ k O ( k ) = (cid:2) D k , O ( k ) (cid:3) . (36)The relation between observables in the velocity andlength gauges can be cast in this language. First, wenote that h ψ k s | O A (cid:18) ∇ i (cid:19) | ψ k ′ s ′ i = h ψ k s | O E (cid:18) ∇ i + e ~ A ( t ) (cid:19) | ψ k ′ s ′ i . (37)Using the form of the Bloch functions, Eq. (10), thisintegral over all space can reduced to one over a unitcell, summed over all of them, giving,(2 π ) δ ( k − k ′ ) h u k s | O A ( k ) | u k s ′ i (38)= (2 π ) δ ( k − k ′ ) h u k s | O E (cid:16) k + e ~ A ( t ) (cid:17) | u k s ′ i (39)or, O A ( k ) = O E (cid:16) k + e ~ A ( t ) (cid:17) . (40)The relationship between operators in the two descrip-tions can also be expressed in terms of these objectsdefined in the same k -point. To do so, we expand theRHS of that equation in powers of A ( t ). It follows fromEq. (36), that O A ( k , t ) = O E (cid:16) k + e ~ A ( t ) , t (cid:17) = + ∞ X n =0 n ! (cid:16) e ~ (cid:17) n A α ( t ) . . . A α n ( t ) (cid:2) D α k , (cid:2) ..., (cid:2) D α n k , O E ( k , t ) (cid:3) ... (cid:3)(cid:3) . (41)where a sum over the repeated cartesian indexes α j is left implied. To conclude this brief account of the use ofthe covariant derivative, we point out that the canonicalcommutation relation, (cid:2) ˆ r α , ˆ p β (cid:3) = i ~ δ αβ ˆ1 , (42)is expressed in the Bloch basis as (cid:2) D α k , V β ( k ) (cid:3) ss ′ = ~ m e δ αβ δ ss ′ , (43)since ˆ r α = iD α and ˆ p α = m e V α . Eq. (43) can alsobe explicitly derived from the form of the D k and v k matrices.We can now write the general many-body Hamiltoni-ans, (2) and (3), in the Bloch description, H E ( t ) = X ss ′ Z d d k (2 π ) d c † k s [ δ ss ′ ǫ k s + ie E ( t ) · D k ss ′ ] c k s ′ , (44) H A ( t ) = X ss ′ Z d d k (2 π ) d c † k s (cid:20) δ ss ′ (cid:18) ǫ k s + e A ( t )2 m e (cid:19) + e A ( t ) · v k ss ′ ] c k s ′ . (45)as well as their respective current operators, J E ( t ) = − e X ss ′ Z d d k (2 π ) d c † k s v k ss ′ c k s ′ , (46) J A ( t ) = − e X ss ′ Z d d k (2 π ) d c † k s (cid:2) v k ss ′ + δ ss ′ em e A ( t ) (cid:3) c k s ′ . (47)The expression of the current in the velocity gauge is aconsequence of Eqs. (40) and (43). Any component ofvelocity in operator form satisfies the general relationbetween k diagonal observables in both gauges, V α,A ( k ) = V α,E (cid:16) k + e ~ A (cid:17) . An expansion of the right hand side in powers of A produces only two terms of orders A (0) and A (1) . Theremaining terms involve two or more derivatives with re-spect to k , which, following Eq. (36), can be expressedin terms of commutators of the covariant derivative and (cid:2) D β k , V α,E ( k ) (cid:3) . As the latter is proportional to the iden-tity operator, Eq. (43), the commutators are exactlyzero and the higher order terms vanish. V α,A ( k ) = V α,E ( k ) + e ~ A β (cid:2) D β k , V α,E ( k ) (cid:3) (48)= V α,E ( k ) + em A α ˆ1 . (49)In matrix form this is v α,A k ss ′ = v α k ss ′ + δ ss ′ em e A α ( t ) , where v k ss ′ is the velocity matrix in the length gauge. III. THE REDUCED DENSITY MATRIX
The reduced density matrix (RDM), is a matrix inband space, defined by the average of momentum con-serving inter band transitions ( s = s ′ ) and the intraband ( s = s ′ ) transitions [20]. ρ α k ss ′ ( t ) = h c † k s ′ c k s i α . (50)The superscript α now denotes the gauge in which thisobject is computed, α = A, E . The average of any op-erator ˆ A is the trace h ˆ A i α = tr (cid:2) ρ α ˆ A (cid:3) , (51)where, ρ α , the full many-body density matrix is ρ α = X n p n | ψ αn i h ψ αn | , (52)and (cid:8) | ψ αn i (cid:9) is a complete set of state vectors. In theSchr¨odinger picture, the time evolution of ρ α ( t ) is gov-erned by the time-evolution of the state vectors, | ψ n ( t ) i , i ~ ∂ | ψ αn ( t ) i ∂t = H α ( t ) | ψ αn ( t ) i . (53)The equation of motion of the RDM takes the form i ~ ∂ρ α k ss ′ ( t ) ∂t = tr (cid:20) i ~ ∂ρ α ( t ) ∂t c † k s ′ c k s (cid:21) = tr h(cid:2) H α ( t ) , ρ α ( t ) (cid:3) c † k s ′ c k s i = (cid:10)(cid:2) c † k s ′ c k s , H α ( t ) (cid:3)(cid:11) α . (54)The gauge freedom to express the uniform electric fieldis thereby carried into the calculation of the system’sdynamics, which can be described either in terms of ρ E k ss ′ ( t ) or ρ A k ss ′ ( t ).Computing the commutators on the right hand sideof Eq. (54), using the Hamiltonians, Eqs. (44) and (45),we obtain closed equations of motion for the RDM [21], (cid:2) i ~ ∂∂t − ǫ k ss ′ (cid:3) ρ E k ss ′ ( t ) = ie E ( t ) · (cid:2) D , ρ E ( t ) (cid:3) k ss ′ , (55) (cid:2) i ~ ∂∂t − ǫ k ss ′ (cid:3) ρ A k ss ′ ( t ) = e A ( t ) · (cid:2) v , ρ A ( t ) (cid:3) k ss ′ . (56)Here we have defined ǫ k ss ′ := ǫ k s − ǫ k s ′ . The equa-tion for the scalar potential RDM is found in references[4, 6, 12], but not, as here, cast in terms of the covariantderivative. The presence of the derivative with respectto k in Eq. (55) couples the response at different valuesof k , whereas, its counterpart for the vector potentialgauge, Eq. (56), is completely decoupled in crystal mo-mentum, k , and can thus be solved independently foreach point of the FBZ. Averages of single particle ob-servables, diagonal in momentum space, such as the cur-rents, (46) and (47), can be obtained from the RDM’sas traces over band space (cid:10) J E ( t ) (cid:11) = − e Z d d k (2 π ) d Tr h V E ( k ) ρ E ( k , t ) i , (57) (cid:10) J A ( t ) (cid:11) = − e Z d d k (2 π ) d Tr h(cid:16) V E ( k ) + em e A ( t )ˆ1 (cid:17) ρ A ( k , t ) i , (58)for ˆ1 the identity in band space. Given that these twoalternative formulations are related by a unitary trans-formation, Eq. (57) and Eq. (58) have to yield the sameresults, although this is far from obvious at this point.To relate these two objects, we must first establish therelation between RDMs.The full many-body density matrix ρ A , ρ A = X n p n (cid:12)(cid:12) ψ An ( t ) (cid:11) (cid:10) ψ An ( t ) (cid:12)(cid:12) , (59)can be expressed in the state vectors of its counter-part description by means of a unitary transformation,Eq. (5), ρ A = U † ( t ) ρ E U ( t ) . (60)The vector potential RDM is therefore expressible asaverages with ρ E , of suitably modified operators ρ A k ss ′ ( t ) = tr h ρ A c † k s ′ c k s i = tr h U † ( t ) ρ E U ( t ) c † k s ′ c k s i = tr h ρ E U ( t ) c † k s ′ c k s U † ( t ) i = tr h ρ E ˜ c † k s ′ ˜ c k s i . (61) The new creation operator ˜ c † k s ′ is obtained from c † k s ′ bythe same unitary transformation˜ c † k s ′ = U ( t ) c † k s ′ U † ( t )= Z d d r Φ k s ′ ( r ) Ψ † ( r ) (62)and creates an electron with wave functionΦ k s ′ ( r ) := e ie r · A ( t ) / ~ ψ k s ′ ( r ) = e i r · ( k + e A ( t ) / ~ ) u k s ′ ( r ) . (63)While this is clearly a Bloch state with wave vector q = k + e A ( t ) / ~ , it is not ψ k + e A ( t ) / ~ s ′ , but can be expandedas a linear combination of the shifted Bloch states, | Φ k s ′ i = X r ′ Z d d q (2 π ) d (cid:12)(cid:12) ψ q r ′ (cid:11)(cid:10) ψ q r ′ (cid:12)(cid:12) Φ k s ′ (cid:11) = X r ′ (cid:12)(cid:12) ψ k + e A ( t ) / ~ r ′ (cid:11)(cid:10) u k + e A ( t ) / ~ r ′ (cid:12)(cid:12) u k s ′ (cid:11) . (64)This allows us to relate ˜ c † k s ′ and ˜ c k s ′ to the originalcreation and destruction operators˜ c † k s ′ = X r ′ (cid:10) u k + e A ( t ) / ~ r ′ (cid:12)(cid:12) u k s ′ (cid:11) c † k + e A ( t ) / ~ r ′ , (65)and express the RHS of Eq. (61) in terms of the shiftedBloch operators. Sincetr (cid:2) ρ E c † k + e A ( t ) / ~ r ′ c k + e A ( t ) / ~ r (cid:3) = h c † k + e A ( t ) / ~ r ′ c k + e A ( t ) / ~ r i E = ρ E k + e A ( t ) / ~ rr ′ ( t ) , (66)using this and Eq. (65) in Eq. (61) we obtain the relationbetween ρ Ass ′ ( k , t ) and ρ Ess ′ ( k , t ) as ρ A k ss ′ ( t ) = X rr ′ h u k s | u k + e A ( t ) / ~ r ′ i ρ E k + e A ( t ) / ~ rr ′ ( t ) h u k + e A ( t ) / ~ r | u k s ′ i . (67)Following our definition of operators in band space,Eq. (29), this equality can be cast in operator form[22],ˆ ρ A ( k , t ) = ˆ ρ E (cid:16) k + e ~ A ( t ) , t (cid:17) . (68)The simplicity of this relation between the density ma-trices in the velocity and length gauges, which amountsto a simple shift in the value of k , is only true in theoperator representation; as can be seen in Eq. (67), thematrix elements of these two operators between statesof a single basis {| u k s i , s = 0 , , . . . } with the same k ,do not satisfy this simple relation.We can now show that the expectation values of ob-servables in the two descriptions are exactly the same.Consider, hO A ( t ) i = Z d d k (2 π ) d Tr (cid:2) O A ( k ) ρ A ( k , t ) (cid:3) . (69) Both the RDM and the observable operators are equal totheir scalar potential counterparts when their argumentis adequately translated, Eqs. (40) and (68), hO A ( t ) i = Z d d k (2 π ) d Tr h O E (cid:16) k + e ~ A (cid:17) ρ E (cid:16) k + e ~ A (cid:17)i . (70)This shifts the FBZ by a constant, which is irrelevantas the integrand is periodic. This means that, hO A ( t ) i = Z d d k (2 π ) d Tr (cid:2) O E ( k ) ρ E ( k , t ) (cid:3) = hO E ( t ) i . (71)The sum rules in reference [12] can be traced back tothis equivalence (see Appendix A). IV. SOLUTIONS TO THE RDM EQUATIONSOF MOTION
A system’s current response (in either scalar or vec-tor potential formalism) is obtained from the solutionto its respective RDM equation of motion. To writethese solutions explicitly, we must employ the generalprocedure from nonlinear physics: break the RDM intocontributions of different powers on the external field, ρ = ρ (0) + ρ (1) + ( ... ) and then proceed to iterativelysolve the equations of motion for each order. For ρ ( n ) ,these read as (cid:20) i ~ ∂∂t − ǫ k ss ′ (cid:21) ρ ( n ) ,E k ss ′ ( t ) = ie E ( t ) · (cid:2) D , ρ ( n − ,E ( t ) (cid:3) k ss ′ , (72) (cid:20) i ~ ∂∂t − ǫ k ss ′ (cid:21) ρ ( n ) ,A k ss ′ ( t ) = e A ( t ) · (cid:2) v , ρ ( n − ,A ( t ) (cid:3) k ss ′ . (73)In expressing the time-dependent objects in termsof their Fourier decomposition, we assume adiabaticswitching of the perturbation ( η → + ), ρ ( n ) ,α k ss ′ ( t ) = Z dω π e − i ( ω + iη ) t ρ ( n ) ,α k ss ′ ( ω ) , (74) E ( t ) = Z dω π e − i ( ω + iη ) t E ( ω ) , (75) A ( t ) = Z dω π e − i ( ω + iη ) t A ( ω ) . (76)The time derivative in the equations of motion is re-placed by a frequency factor that is collected into an energy denominator, d k ss ′ ( ω ) := 1 ~ ω − ǫ k ss ′ . (77)From this point on, the frequency argument of theseenergy denominators is understood to have an infinites-imal imaginary part. Using the Hadamard product oftwo matrices, ( A ◦ B ) ss ′ := A ss ′ B ss ′ , (78)we write recursion relations for the RDM solutions, ρ ( n ) ,E k ss ′ ( ω ) = ie Z dω π E α ( ω ) ( d ( ω ) ◦ (cid:2) D α , ρ ( n − ,E ( ω − ω ) (cid:3) ) k ss ′ , (79) ρ ( n ) ,A k ss ′ ( ω ) = e Z dω π A α ( ω ) ( d ( ω ) ◦ (cid:2) v α , ρ ( n − ,A ( ω − ω ) (cid:3) ) k ss ′ . (80)The successive application of these expressions, bringsthe n -th order solution, ρ ( n ) k ss ′ ( ω ), to the form of nestedcommutators of the zeroth order one, ρ (0) k ss ′ , which is theFermi-Dirac distribution function times the unit matrixin band space, ρ (0) k ss ′ = f k s δ ss ′ . With one last bit ofnotation, ω [ m ] := m X i =1 ω i , (81)we write ρ ( n ) k ss ′ in each formalism as ρ ( n ) ,E k ss ′ ( ω ) = ( ie ) n " n − Y i =1 Z dω i π E α i ( ω i ) E α n ( ω − ω [ n − ) (cid:16) d ( ω ) ◦ (cid:2) D α , d ( ω − ω ) ◦ (cid:2) ..., d ( ω − ω [ n − ) ◦ (cid:2) D α n , ρ (0) (cid:3) ... (cid:3)(cid:3)(cid:17) k ss ′ , (82) ρ ( n ) ,A k ss ′ ( ω ) = e n " n − Y i =1 Z dω i π A α i ( ω i ) A α n ( ω − ω [ n − ) (cid:16) d ( ω ) ◦ (cid:2) v α , d ( ω − ω ) ◦ (cid:2) ..., d ( ω − ω [ n − ) ◦ (cid:2) v α n , ρ (0) (cid:3) ... (cid:3)(cid:3)(cid:17) k ss ′ . (83)The use of the nested commutator and the factor d ( ω )allow for a compact form of these solutions (albeit hid-ing their considerable complexity). These solutions areentirely written in terms of three objects (and theirderivatives, in the case of ρ E ): the band energies, ǫ k s ; the Berry connection, ξ k ss ′ ; and the structure offilled/empty bands, f k s . In principle, by determiningthese, one determines a system’s nonlinear response, provided that one can compute the FBZ integrals inEqs.(57) and (58). In section VI, we use these expres-sions to discuss whether it is possible to truncate thesums over bands, implicit in these matrix products, toa reduced set. But before, we study the first two non-trivial orders of the monolayer graphene response to auniform electric field, to illustrate that these reproducethe results of references [4, 6, 8]. V. LINEAR AND THIRD ORDER RESPONSEIN THE MONOLAYER GRAPHENE
Monolayer graphene is (usually) described by a tight-binding model that considers only nearest neighbourhopping[23] and for which the expressions for the twobands, ǫ k s , and the periodic functions, u k s , can be ex-plicitly computed. This allows one to derive some usefulproperties.First, a double band-index sum, such as the ones in h J E i and h J A i , over an antisymmetric object θ ss ′ = − θ s ′ s , can be reduced to a single sum, where ¯ s reads asthe band opposite to s , X s ′ s θ ss ′ → X s θ s ¯ s , (84)Second, the Berry connection for the monolayer is de-scribed by an intra band and an inter band term. Withadequate choice of gauge (Eq. 24) the following proper-ties can be obtained [4, 6, 8], ξ α k ss = ξ α k ¯ s ¯ s , (85) ξ α k s ¯ s = ξ α k ¯ ss . (86)In addition, we can choose these to be even under k →− k . Linear order response
The first order term of the RDM is ρ (1) k ss ′ ( ω ) = ie E α ( ω ) (cid:16) d ( ω ) ◦ (cid:2) D α , ρ (0) (cid:3)(cid:17) k ss ′ , (87)= ie E α ( ω ) n δ ss ′ ~ ω ∇ α k f k s − i ξ α k ss ′ ( f k s ′ − f k s ) ~ ω − ǫ k ss ′ o . (88)Writing the current, Eq. (57), in terms of its Fouriercomponents, we obtain two contributions [4, 6, 8], h J (1) ,β ( ω ) i = E α ( ω ) Z d d k (2 π ) d (cid:2) Π (1) ,βα k , intra ( ω ) + Π (1) ,βα k , inter ( ω ) (cid:3) . (89)The intra band contribution is a generalized Drudeterm,Π (1) ,βα k , intra ( ω ) = i e ~ ω X s (cid:0) ∇ β k ǫ k s (cid:1)(cid:0) ∇ α k ǫ k s (cid:1) (cid:18) − ∂f k s ∂ǫ k s (cid:19) . (90)and the inter band contribution, Π αβ, (1)inter ( k , ω ), whichinvolves the Berry connection,Π (1) ,βα k , inter ( ω ) = − e X s v β k ,s ¯ s ξ α k , ¯ ss f k s − f k ¯ s ~ ω − ǫ k ¯ ss . (91) Third order response
Because graphene has inversion symmetry, its secondorder response is zero (see Appendix B). The first non-linear contribution to the current is thus the third orderone, which is to be computed here. From the generalRDM solution, Eq. (82), we obtain for n = 3, ρ (3) k ss ′ ( ω ) = ( ie ) Z dω π Z dω π E α ( ω ) E α ( ω ) E α ( ω − ω [2] ) (cid:16) d ( ω ) ◦ (cid:2) D α , d ( ω − ω ) ◦ (cid:2) D α , d ( ω − ω [2] ) ◦ (cid:2) D α , ρ (0) (cid:3)(cid:3)(cid:3)(cid:17) k ss ′ . (92)Expanding this by means of Eq. (34), producesa plethora of terms which we organize, followingMikhailov, by the number of intra band (derivatives)and inter band (Berry connections) factors. These canbe manipulated independently and give contributions to the current, that can be labeled by i = 1 , ,
3, the num-ber of intra band factors.The field factors and the FBZ integrations are thesame in all these contributions, and so h J (3) ,β ( ω ) i = e Z dω π Z dω π E α ( ω ) E α ( ω ) E α ( ω − ω [2] ) Z d d k (2 π ) d X i =0 Π (3) ,βα α α k ,i ( ω, ω , ω ) . (93)The Π contribution is the single term in Eq. (92) with three k -derivatives, while Π gathers three terms withtwo k -derivatives, to which we can apply the propertyEq. (84). Contributions Π and Π require additional manipulations, presented in Appendix C. Consistentlywith our notation ǫ k s ¯ s := ǫ k s − ǫ k ¯ s , we abbreviate thefactor f k s − f k ¯ s as f k s ¯ s . The expressions for the Π i are:Π (3) ,βα α α = i ~ ω ~ ( ω − ω ) 1 ~ ( ω − ω [2] ) X s v β k ss ∇ α k ∇ α k ∇ α k f k s , (94)Π (3) ,βα α α = X s v β k ¯ ss ~ ω − ǫ k s ¯ s (cid:26) ~ ( ω − ω ) 1 ~ ( ω − ω [2] ) ξ α k s ¯ s ∇ α k ∇ α k f k ¯ ss + 1 ~ ( ω − ω [2] ) ∇ α k (cid:18) ξ α k s ¯ s ∇ α k f k ¯ ss ~ ( ω − ω ) − ǫ k s ¯ s (cid:19) (95)+ ∇ α k (cid:18) ~ ( ω − ω ) − ǫ k s ¯ s ∇ α k (cid:18) ξ α k s ¯ s f k ¯ ss ~ ( ω − ω [2] ) − ǫ k s ¯ s (cid:19)(cid:19)(cid:27) , (96)Π (3) ,βα α α = i ~ ω X s v β k ss (cid:26) ~ ( ω − ω ) ∇ α k (cid:18) ξ α k s ¯ s ξ α k ¯ ss f k s ¯ s (cid:18) ~ ( ω − ω [2] ) − ǫ k ¯ ss + 1 ~ ( ω − ω [2] ) − ǫ k ¯ ss (cid:19)(cid:19) + ξ α k s ¯ s ~ ( ω − ω ) − ǫ k ¯ ss ∇ α k (cid:18) ξ α k ¯ ss f k s ¯ s (cid:18) ~ ( ω − ω [2] ) − ǫ k ¯ ss + 1 ~ ( ω − ω [2] ) − ǫ k s ¯ s (cid:19)(cid:19) + 1 ~ ( ω − ω [2] ) ξ α k s ¯ s ξ α k ¯ ss (cid:18) ~ ( ω − ω ) − ǫ k ¯ ss + 1 ~ ( ω − ω ) − ǫ k s ¯ s (cid:19) ∇ α k f k s ¯ s (cid:27) , (97)Π (3) ,βα α α = 2 ~ ( ω − ω ) X s ~ ω − ǫ k s ¯ s v β k ¯ ss ξ α k s ¯ s ξ α k ¯ ss ξ α k s ¯ s f k ¯ ss (cid:18) ~ ( ω − ω [2] ) − ǫ k s ¯ s + 1 ~ ( ω − ω [2] ) − ǫ k ¯ ss (cid:19) . (98)Apart from differences in the Cartesian and frequencyindexes, which can always be relabeled, these are theexpressions found in [4, 6, 8], in the limit where Γ ( i ) ,Γ ( e ) , the phenomenological scattering rates, are set tozero. VI. EFFECTIVE HAMILTONIANS
As we have seen in the previous sections, the currentresponse in the Schr¨odinger problem can be computedby two different but equivalent procedures. Althoughthis is conceptually important, one is ultimately inter-ested in computing the current in materials describedby effective Hamiltonians, that account only for a finitenumber of bands.We will now show that, by truncating the band space,the scalar and vector potential currents are no longerthe same; the latter contains relevant contributions frombands that are left out of the effective Hamiltonian.In Fig. (1), we present a conceptual picture of a spec-trum which has a cluster of bands close to the Fermilevel, which we deem relevant, well separated in en-ergy by the bands below (filled) and above (empty). Wedenote the energy scale in the subspace E of relevantbands by δ and the energy separation to other bands by ε ∆∆δ Figure 1. A conceptual picture of an effective Hamiltonianthat describes the bands in subspace E . The Fermi levellies somewhere in that subspace. For bands inside E , theenergy difference is of order δ , while the energy differencebetween bands in different subspaces is of order ∆ ≫ δ . ∆ ≫ δ ; we assume the frequency of the external field, ω , to be of the order, ω ∼ δ/ ~ . The question we wish toanswer is whether the bands outside E can be ignoredin the calculation of the current.Energy denominators, d ( ω ), involving bands inside E , d k ss ′ ( ω ) = 1 ~ ω − ǫ k ss ′ , ǫ k ss ′ ∼ δ, E and those outside d k ss ′ ( ω ) = 1 ~ ω − ǫ k ss ′ , ǫ k ss ′ ∼ ∆ . In the n -th order contribution to the current in thescalar potential gauge, each term is a trace of a product n + 1 matrices in band space. X s ′ s v β k s ′ s d k ss ′ ( ω ) (cid:2) D α , d ( ω − ω ) ◦ (cid:2) ... (cid:2) D α n , ρ (0) (cid:3) ... (cid:3)(cid:3) k ss ′ . (99)One such term is the fully intra band one, in which wepick only the diagonal part of each covariant derivativeoperator (see an example in Eq. (94)). This term is nonzero only for the band that contains the Fermi level; ithas no contributions from bands outside E . All inter band contributions contain at least one dif-ference of occupation factors, thus allowing us to discardany terms that involve only filled, or only empty bands.What remains are contributions of three types: (a)terms that involve transition between filled and emptybands, both outside E ; (b) terms that involve transitionbetween bands in E and bands outside E ; (c) inter bandterms among the bands in E .In Eq. (99), if s ′ belongs to a filled band outside E ,and s to a empty band, also outside E (or vice-versa), ǫ k ss ′ ∼ ∆ ≫ ~ ω and v β k s ′ s d k ss ′ ( ω ) ∼ i ~ ǫ k ss ′ ξ β k s ′ s ~ ω − ǫ k ss ′ = − i ~ ξ β k s ′ s ;the matrix that follows this term in Eq. (99), must in-clude at least an energy denominator ∼ O (1 / ∆) becausethe last band index is the same as the first, s ′ . If s ′ and s both refer to a filled band (or an empty one) outside E the inter band terms between filled and empty bandsmust have at least two such energy denominators. Inother words, terms of type (a) have a least and energydenominator of order ∼ O (1 / ∆) . An identical argumentcan be made for transitions between bands inside andand outside E , i.e. for transitions of type (b). We con-clude that, in the limit of ∆ ≫ δ , the dominant termscome from bands in E (terms of type (c)) and trunca-tion to this subspace is a valid approximation.The same argument does not carry to the correspond-ing contribution to the current in the velocity gauge: X s ′ s v β k s ′ s d k ss ′ ( ω ) (cid:2) v α , d ( ω − ω ) ◦ (cid:2) ... (cid:2) v α n , ρ (0) (cid:3) ... (cid:3)(cid:3) k ss ′ . (100)In this case, every energy denominator is associated witha velocity matrix element. In any transition involvingenergies ǫ rr ′ ∼ O (∆), the d rr ′ ( ω ) ∼ O (1 / ∆), as before,but the corresponding velocity matrix element has anoff-diagonal contribution v β k r ′ r ∼ ξ β k rr ′ O (∆), and suchterms give relevant contributions no matter how large∆ is. In other words, bands away from the Fermi surfacecontribute just as much as those in E for the expressionof the current in the velocity gauge. VII. SUMMARY AND CONCLUSIONS
The concept of covariant derivative in k -space hasbeen shown here to be of considerable value in thecalculation of the nonlinear current response. It is avery convenient representation of the position operator[Eq. (23)]; it clarifies the structure of the velocity matrix[Eq. (27)]; it allows a complete parallel development ofthe structure of the reduced density matrices (RDM) inthe length (scalar potential) and velocity (vector poten-tial) gauges [Eqs. (55) and (56)] and it provides compactexpressions for the perturbative solutions of the equa-tions of motion of the RDM [Eqs. (82) and (83)]. It alsoallowed us to see how the equivalence between objectsin the two gauges breaks down when the band space istruncated.Furthermore that the approximation is legitimate inthe length gauge, but fails in the velocity one is madeclear by considering the truncation of the covariantderivative and the velocity operators to a restricted setof bands. The commutator of these two quantities,Eq. (43), which is constant for the case of the infinitebands of the Schr¨odinger Hamiltonian, no longer holdsfor a truncated subset of these. In fact, the statementthat the commutator is constant is equivalent to theHamiltonian being linear or quadratic in k (the formeris the case of the Dirac Hamiltonian, where the com-mutator gives zero). As a result, Eq. (48) is no longervalid for a general effective Hamiltonian and there is noequality between currents in the two gauges, at least asthey are written in Eqs. (46) and (47). In order to usethe velocity gauge in actual calculations, one must startfrom the beginning with the effective Hamiltonian andperform the minimum coupling then. Naturally, thismeans modifying both the equation of motion of theRDM, Eq. (56), and the current operator in the velocitygauge. This will be the subject of a future paper. VIII. ACKNOWLEDGMENTS
The work of G.B.V and D.J.P is supported by Fun-da¸c˜ao para a Ciˆencia e Tecnologia (FCT) under thegrants PD/BI/129220/2017 and PD/BD/135019/2017respectively. N.M.R.P. acknowledges funding from theEuropean Commission within the project “Graphene-Driven Revolutions in ICT and Beyond” (ref. no.696656) and the Portuguese Foundation for Science andTechnology (FCT) in the framework of the Strategic Fi-nancing UID/FIS/04650/2013.
Appendix A: The Aversa and Sipe Sum rules
In reference [12], Aversa and Sipe sketch how prob-lems might arise in the velocity gauge as it involves ad-ditional terms (with respect to the length gauge) thatare zero only when they are treated exactly. They showthis explicitly in the linear response.1In our formalism this can be seen in full generality,starting from Eqs. (41) and (70). Expanding in powersof A ( t ), we can write hO A ( k , t ) i in terms of hO E ( k , t ) i and contributions that depend explicitly on A ( t ), thatis, the n ≥ hO A ( k , t ) i = hO E ( k , t ) i + e ~ A α ( t ) Z d d k (2 π ) d Tr (cid:2) D α , O E ( k , t ) ρ E ( k , t ) (cid:3) +( ... ) , (A1)The equivalence of the two gauges requires that every-thing other than hO E ( k , t ) i in the right hand side tobe zero. Each contribution contains a factor which, forarbitrary order n , reads as Z d d k (2 π ) d Tr (cid:2) D α , G( k , t ) (cid:3) , (A2)for G( k , t ) is some matrix in band space,G( k , t ) := (cid:2) D α , ..., (cid:2) D α n , (cid:2) O E ( k , t ) ρ E ( k , t ) (cid:3) ... (cid:3)(cid:3) , The two terms in the expression (A2) amount to X s Z d d k (2 π ) d ∇ k G ss ( k , t ) − i Z d d k (2 π ) d Tr [ ξ α , G( k , t)] = 0 . The first term integrates to zero due to the periodicityof G( k , t ) in reciprocal space. As for the second one, thetrace of any commutator of two matrices is zero, whichfollows from the cyclic invariance of the trace.The cancellation of all terms in the sum on the RHS ofEq. (A1) constitute the sum rules referred to by Aversaand Sipe. They are an order by order formulation of theresult that in an integral over the FBZ of a function withthe period of the reciprocal lattice, the argument of theintegrand can be shifted without changing the integral.Nevertheless, it is clear from this formulation thatthese extra terms exist only if, starting from the velocitygauge, we try to reduce our expressions to the ones inthe length gauge. Appendix B: Current response in acentrosymmetric material
For a centrosymmetric crystal, the spatial inversionoperator P commutes with the crystal Hamiltonian, H ,[ P , H ] = 0 . (B1)and the solutions to the k -dependent Hamiltonian, u k s ( − r ) and u − k s ( r ), are related by a phase phase fac-tor, u k s ( − r ) = e iµ k s u − k s ( r ) . (B2)If we take r → − r in the integral that defines the Berryconnection, Eq. (21), we can determine how it behaveswhen we exchange the sign of the crystal momentum, ξ k ss ′ = iv C Z uc d d r u ∗ k s ( − r ) ∇ k u k s ′ ( − r )= − e i ( µ k s ′ − µ k s ) (cid:2) ξ − k ss ′ + ∇ k µ k s ′ δ ss ′ (cid:3) , (B3) and, ξ − k ss ′ = − e i ( µ k s − µ k s ’ ) [ ξ k ss ′ + ∇ k µ k s ′ δ ss ′ ] . (B4)This determines the transformation law for the covariantderivative: D − k ss ′ = − e i ( µ k s − µ k s ’ ) [ D k ss ′ + iδ ss ′ ∇ k µ k s ′ ] . (B5)We can now determine how the k -dependent factor of ρ (1) , Π (1) ,α k ss ′ ( ω ) = d k ss ′ ( ω ) (cid:2) D α , ρ (0) (cid:3) k ss ′ , (B6)transforms upon k → − k , by recalling that f k s and d k ss ′ ( ω ) are even, since they depend only on the bandenergies,Π (1) ,α − k ss ′ ( ω ) = d − k ss ′ ( ω ) (cid:2) D α , ρ (0) (cid:3) − k ss ′ = d k ss ′ ( ω ) (cid:0) D α − k ss ′ f k s ′ − f k s D α − k ss ′ (cid:1) = − e i ( µ k s − µ k s ’ ) d k ss ′ ( ω ) (cid:2) D α , ρ (0) (cid:3) k ss ′ = − e i ( µ k s − µ k s ’ ) Π (1) ,α k ss ′ ( ω ) . (B7)This can be extended to the higher order contributionsof the RDM, in particular to the second order one, easilyextracted from Eq. (82),Π (2) ,α α k ss ′ ( ω, ω ) = d k ss ′ ( ω ) (cid:2) D α , Π (1) ,α ( ω − ω ) (cid:3) k ss ′ . (B8)It follows from Eqs. (B5) and (B7) thatΠ (2) ,α α − k ss ′ ( ω, ω ) = e i ( µ k s − µ k s ’ ) Π (2) ,α α k ss ′ ( ω, ω ) . (B9)This object picks up the same k -space phase factors,but unlike its first order counterpart, the sign does notchange. Combining this with the transformation law forthe velocity matrix element, v k s ′ s → v − k s ′ s = − e i ( µ k s ′ − µ k s ) v k s ′ s . (B10)we see that the integrand in the FBZ integral, v β k s ′ s Π (2) ,α α k ss ′ ( ω, ω ) is an odd function of k , the sec-ond order current vanishes.This argument also carries for an arbitrary order n .The k -parity of v β k s ′ s Π ( n ) ,α ( ... ) α n k ss ′ is determined by itsnumber of covariant derivatives in Π ( n ) ,α ( ... ) α n k ss ′ . For n even, the integrand is odd under k → − k , so even ordercontributions to the current vanish in a centrosymmetricmaterial. Appendix C: Deriving the expressions: (97)-(98)
Consider Π , the collection of terms with only oneintra band factor, where we have used the two-bandcharacter of the monolayer graphene, Eq. (84),2Π (3) ,βα α α = i X r s ′ s v β k s ′ s (cid:26) ~ ω − ǫ k ss ′ ∇ α k (cid:18) ~ ( ω − ω ) − ǫ k ss ′ (cid:18) δ r ¯ s ′ ξ α k s ¯ s ′ ξ α k ¯ s ′ s ′ f k s ′ ¯ s ′ ~ ( ω − ω [2] ) − ǫ k ¯ s ′ s ′ − δ r ¯ s ξ α k s ¯ s ξ α k ¯ ss ′ f k ¯ ss ~ ( ω − ω [2] ) − ǫ k s ¯ s (cid:19)(cid:19) + 1 ~ ω − ǫ k ss ′ (cid:18) δ r ¯ s ′ ξ α k s ¯ s ′ ~ ( ω − ω ) − ǫ k ¯ s ′ s ′ ∇ α k (cid:18) ξ α k ¯ s ′ s ′ f k s ′ ¯ s ′ ~ ( ω − ω [2] ) − ǫ k ¯ s ′ s ′ (cid:19) − δ r ¯ s ξ α k ¯ ss ′ ~ ( ω − ω ) − ǫ k s ¯ s ∇ α k (cid:18) ξ α k s ¯ s f ¯ ss ~ ( ω − ω [2] ) − ǫ k s ¯ s (cid:19)(cid:19) + 1 ~ ω − ǫ k ss ′ (cid:18) δ r ¯ s ′ ξ α k s ¯ s ′ ξ α k ¯ s ′ s ′ ~ ( ω − ω ) − ǫ k ¯ s ′ s ′ ∇ α k f k s ′ ¯ s ′ ~ ( ω − ω [2] ) − δ r ¯ s ξ α k s ¯ s ξ α k ¯ ss ′ ~ ( ω − ω ) − ǫ k s ¯ s ∇ α k f k ¯ ss ~ ( ω − ω [2] ) (cid:19)(cid:27) . (C1)For s ′ = ¯ s , the two terms in each line cancel out byapplication of the Berry connection properties, Eq. (85) and (86). This fixes s ′ = s , and the Π contributionreads asΠ (3) ,βα α α = i ~ ω X s v β k ss n ~ ( ω − ω ) ∇ α k (cid:18) ξ α k s ¯ s ξ α k ¯ ss f k s ¯ s (cid:18) ~ ( ω − ω [2] ) − ǫ k ¯ ss + 1 ~ ( ω − ω [2] ) − ǫ k s ¯ s (cid:19)(cid:19) + ξ α k s ¯ s ~ ( ω − ω ) − ǫ k ¯ ss ∇ α k (cid:18) ξ α k ¯ ss f k s ¯ s (cid:18) ~ ( ω − ω [2] ) − ǫ k ¯ ss + 1 ~ ( ω − ω [2] ) − ǫ k s ¯ s (cid:19)(cid:19) (C2)+ 1 ~ ( ω − ω [2] ) ξ α k s ¯ s ξ α k ¯ ss (cid:18) ~ ( ω − ω ) − ǫ k ¯ ss + 1 ~ ( ω − ω ) − ǫ k s ¯ s (cid:19) ∇ α k f k s ¯ s o . (C3)The Π portion of the current for the two band material reduces to,Π (3) ,βα α α = X r ′ r s ′ s v β k s ′ s ~ ω − ǫ k ss ′ (cid:26) ~ ( ω − ω ) − ǫ k rs ′ δ r ′ ¯ s ′ ξ α k sr ξ α k r ¯ s ′ ξ α k ¯ s ′ s ′ f k s ′ ¯ s ′ ~ (cid:0) ω − ω [2] (cid:1) − ǫ k ¯ s ′ s ′ − δ r ′ ¯ r ξ α k sr ξ α k r ¯ r ξ α k ¯ rs ′ f k ¯ rr ~ (cid:0) ω − ω [2] (cid:1) − ǫ k r ¯ r ! − ~ ( ω − ω ) − ǫ k sr δ r ′ ¯ r ξ α k s ¯ r ξ α k ¯ rr ξ α rs ′ f k r ¯ r ~ (cid:0) ω − ω [2] (cid:1) − ǫ k ¯ rr − δ r ′ ¯ s ξ α k s ¯ s ξ α k ¯ sr ξ α rs ′ f k ¯ ss ~ (cid:0) ω − ω [2] (cid:1) − ǫ k s ¯ s !o . (C4)Since the Berry connection is even under k → − k , thevelocity matrix element in band space, Eq. (28), is writ-ten as the sum of two contributions of opposite parity:an odd intra band term and an even inter band term. The integration over the FBZ carries cancels all oddterms, and so the intra band part of v k ss ′ can be ig-nored: this fixes s ′ = ¯ s , andΠ (3) ,βα α α = X r s v β k ¯ ss ~ ω − ǫ k s ¯ s (cid:26) ~ ( ω − ω ) − ǫ k r ¯ s ξ α k sr ξ α k rs ξ α k s ¯ s f k ¯ ss ~ (cid:0) ω − ω [2] (cid:1) − ǫ k s ¯ s − ξ α k sr ξ α k r ¯ r ξ α k ¯ r ¯ s f k ¯ rr ~ (cid:0) ω − ω [2] (cid:1) − ǫ k r ¯ r ! − ~ ( ω − ω ) − ǫ k sr ξ α k s ¯ r ξ α k ¯ rr ξ α k r ¯ s f k r ¯ r ~ ω − ǫ k ¯ rr − ξ α k s ¯ s ξ α k ¯ sr ξ α r ¯ s f k ¯ ss ~ (cid:0) ω − ω [2] (cid:1) − ǫ k s ¯ s !(cid:27) . (C5)Setting r = s and using (85) and (86), the first twoterms of this expression cancel out; the last two cancel when r = ¯ s .Finally, Π reduces to,3Π (3) ,βα α α = 2 ~ ( ω − ω ) X s ~ ω − ǫ k s ¯ s v β k ¯ ss ξ α k s ¯ s ξ α k ¯ ss ξ α k s ¯ s f k ¯ ss ( ~ (cid:0) ω − ω [2] (cid:1) − ǫ k s ¯ s + 1 ~ (cid:0) ω − ω [2] (cid:1) − ǫ k ¯ ss ) . (C6) [1] S. A. Mikhailov and K. Ziegler,Journal of Physics: Condensed Matter , 384204 (2008).[2] F. Bonaccorso, Z. Sun, T. Hasan, and A. C. Ferrari,Nature Photonics , 611 (2010).[3] M. M. Glazov and S. D. Ganichev, Physics Reports ,101 (2014).[4] J. L. Cheng, N. Vermeulen, and J. E. Sipe,New Journal of Physics , 053014 (2014).[5] N. M. R. Peres, Y. V. Bludov, J. E. Santos, A.-P. Jauho,and M. I. Vasilevskiy, Phys. Rev. B , 125425 (2014).[6] J. L. Cheng, N. Vermeulen, and J. E. Sipe,Phys. Rev. B , 235320 (2015).[7] T. G. Pedersen, Phys. Rev. B , 235432 (2015).[8] S. A. Mikhailov, Phys. Rev. B , 085403 (2016).[9] D. J. Moss, E. Ghahramani, J. E. Sipe, and H. M. VanDriel, Phys. Rev. B , 15 (1990).[10] J. E. Sipe and E. Ghahramani,Phys. Rev. B , 11705 (1993).[11] V. N. Genkin and M. P. Mednis,Sov. Phys. JETP , 609 (1968).[12] C. Aversa and J. E. Sipe, Phys. Rev. B , 14636 (1995).[13] W. E. Lamb, R. R. Schlicher, and M. O. Scully,Phys. Rev. A , 2763 (1987). [14] W. E. Lamb, Phys. Rev. , 259 (1952).[15] K. Rz ֒ a˙zewski and R. W. Boyd, Journal of Modern Op-tics , 1137 (2004).[16] The charge q has already been replaced by the chargeof an electron, q = − e .[17] N. Ashcroft and N. Mermin, Solid State Physics , HRWInternational Editions (Holt, Rinehart and Winston,1976).[18] E. I. Blount,
Solid State Physics: Advances in Researchand Applications , edited by F. Seitz and D. Turnbull,Vol. 13 (Academic, New York, 1962).[19] L. Pitaevskii and E. Lifshitz,
Statistical Physics, Part2. Vol.9 (Butterworth-Heinemann, 1980).[20] Note the switch in the band indexes of the RDM.[21] In order to simplify notation we will drop the crystalmomentum label from the objects inside the commuta-tors and write it alongisde the band indexes.[22] The hats will be henceforth dropped.[23] A. H. Castro Neto, F. Guinea, N. M. R.Peres, K. S. Novoselov, and A. K. Geim,Rev. Mod. Phys.81