Bulk photovoltaic effects in the presence of a static electric field
BBulk photovoltaic effects in the presence of a static electric field
Benjamin M. Fregoso
Department of Physics, Kent State University, Kent, Ohio 44242, USA
This paper presents a study of dc photocurrents in biased insulators to the third order in theelectric field. We find three photocurrents which are characterized by physical divergences of thethird-order free-electron polarization susceptibility. In the absence of momentum relaxation andsaturation effects, these dc photocurrents grow as t n ( n = 2 , ,
0) with illumination time. Thephotocurrents are dubbed jerk , third-order injection, and third-order shift current, respectively, andare generalizations of the second-order injection and shift currents of the bulk photovoltaic effect.We also revisit the theory of the bulk photovoltaic effect and include Fermi surface contributionswhich are important in metals. Finally, we show that injection, shift, and jerk currents admitsimple physical interpretations in terms of semiclassical wave packet dynamics in electric fields.Experimental signatures and extensions to higher-order susceptibilities are also discussed.
I. INTRODUCTION AND MAIN RESULTS
Electrons in crystals can exhibit fascinating dynamicsin the presence of external electric and magnetic fields. Inmetals, the anomalous Hall effect or the chiral anomalyin Weyl semimetals are two examples. Insulators, de-spite lacking a Fermi surface, can also exhibit nontriv-ial carrier dynamics as in the bulk photovoltaic effects(BPVE).
The BPVE is the generation of a dc pho-tocurrent in homogeneous insulators or semiconductorsthat lack inversion symmetry.In the presence of a static E and an optical field E ,the dc current can be expanded in powers of the electricfields. Schematically we can write J dc = σ (1) dark E + σ (2) dark E + σ (2) bpve E + σ (3) ph E E + · · · . (1)The first and second terms are the linear and quadratic dcconductivity in the absence of illumination. The third isthe BPVE and the fourth is the photoconductivity, i.e.,the intensity-dependent dc conductivity. In insulators,the BPVE is usually the dominant contribution. In thisarticle we first review the theory of the BPVE (includingFermi surface contributions), and then we extend it tostudy the photoconductivity.The peculiar nature of the BPVE was first noticedby (1) its dependence on the intensity of light, (2) itslarge open-circuit photovoltages, and (3) its dependenceon light polarization . (1) indicates that the BPVEis quadratic in the optical field, (2) indicates that theBPVE is an ultrafast effect in which transport occursbefore carriers pretermalize at the bottom of the conduc-tion band (top of the valence band), and (3) indicatesthat the BPVE response tensor is complex and has twocomponents. The real part σ couples to the real electricfields and the imaginary part η to the imaginary electricfields, schematically J (2) dc,pbve = σ | E | + η E × E ∗ . (2) This led to the first successful phenomenological theoryof both components of the BPVE, namely, the injectioncurrent, also called circular photogalvanic effect (CPGE),represented by η and shift current represented by σ .The lack of inversion symmetry could manifest in twodistinct scenarios in the BPVE. In the first scenario, pho-toexcited carriers relax momentum asymmetrically into ± k directions via collisions with other electrons, phononsor impurities. This leads to a polar distribution and anet current . In the second scenario, the origin of theBPVE is the light-matter interactions not the dynamicsof momentum relaxation. In injection current processes,light pumps carriers into velocity-carrying states asym-metrically at ± k points in the Brillouin zone (BZ) leadingto a polar distribution and a net current . Within asimple relaxation time approximation, the steady stateinjection current is proportional to the first power of therelaxation time constant and vanishes for linearly polar-ized light.In shift current processes, inversion symmetry break-ing manifests as a separation of the centers of charge ofthe valence and conduction bands so that charge movescoherently across the unit cell upon carrier photoexcita-tion from valence to conduction band . The shift cur-rent vanishes for circular polarization of light and decaysin the time scale of the quantum coherence of the solid.The BPVE has been extensively studied sincethe 1960s in ferroelectrics mainly in the contextof photovoltaic applications. The injection cur-rent, the shift current, or both have been reportedin many materials, includingGaAs, CdSe,
CdS, quantum wells, RhSi, and Bi GeO . More recently, the BPVE has at-tracted attention for its promise in novel optoelectricapplications; specifically in two-dimensional (2D)ferroelectrics.
Following Sipe and coworkers , BPVE response ten-sors can be derived from the perspective of divergent po-larization susceptibilities. In this approach, the BPVEarises from light-matter interactions and not from mo-mentum relaxation processes; the latter are included phe-nomenologically a posteriori . For not too large electricfields, the insulator’s response to an external electric field a r X i v : . [ c ond - m a t . m e s - h a ll ] A ug TABLE I. Summary of bulk photovoltaic effects (BPVEs) obtained from divergences of free electric polarization susceptibilities.The standard BPVEs are derived from the singularities of χ . Higher-order BPVEs can be classified by their dependence onillumination time in the absence of momentum relaxation and saturation effects, e.g., η , η and η are all injection currentresponses and σ , σ are all shift currents responses. We write susceptibilities as, χ abc...n ( − ω Σ , ω β , ω σ , ... ) where b, c... areCartesian indices, ω β , ω σ , ... are frequency components, and ω Σ = ω β + ω σ + ... frequency sums . [ X, Y ] ( { X, Y } ) indicatecommutation (anticommutation) with respect to b, c indices only. Other conventions are explained in Sec. II. BPVE Symbol Expression Timedepen-dence ∼ t α Origin Ref.Injection η abc πe (cid:126) V (cid:80) nm k f mn ω nm ; a [ r bnm , r cmn ] δ ( ω nm − ω ) 1 χ (0 , ω, − ω ) → ∞ σ abc iπe (cid:126) V (cid:80) nm k f mn { r cnm ; a , r bmn } δ ( ω nm − ω ) 0 42Jerk ι abcd πe (cid:126) V (cid:80) nm k f mn (cid:2) ω nm ; ad r bnm r cmn + ω nm ; a ( r bnm r cmn ) ; d (cid:3) δ ( ω nm − ω ) 2 χ (0 , ω, − ω, → ∞ η Eq. 137 1 presentShift σ Eq. 161 0 presentInjection ν Eq. 7 1 χ (0 , − ω, ω, ω ) → ∞ σ Eq. 8 0 44Snap ς Eq. 174 3 χ (0 , ω, − ω, , → ∞ presentJerk ι η σ a n πi (cid:72) | z | = ρ dz χnzl +1 , z = − iω Σ , ρ → , l = − n, · · · , −
1, Eq. 182 α = n − , ..., χ n ( ω Σ , ω β , ω σ , ... ) → ∞ , ω Σ = ω β + ω σ + · · · → is described perturbatively by susceptibilities χ n as P = P + χ E + χ E + χ E + · · · , (3)where P is the electric polarization in the absence of anexternal electric field, χ is the linear susceptibility,and χ , χ , ... are nonlinear susceptibilities. The electric polarization in insulators is commonlythought to be determined by the off-diagonal elementsof the density matrix because these elements describethe displacement of charge from its equilibrium posi-tion in the presence of an electric field. Intraband pro-cesses, however, have been shown to be important.
Among other things they cure unphysical divergences insusceptibilities in the dc limit by incorporating the factthat the intraband motion of Bloch electrons cannot ac-celerate indefinitely in insulators . Importantly, whenintraband and interband processes are taken into accounton an equal footing divergent susceptibilities representreal photocurrents.Consider, for example, the dc divergences of χ . Ifwe denote the amplitude of the electric field by E b = (cid:80) β E bβ e − iω β t , the polarization to second order P a (2) = (cid:88) bβcσ χ abc ( − ω Σ , ω β , ω σ ) E bβ E cσ e − iω Σ t , (4)oscillates with frequency ω Σ = ω β + ω σ in the long-timelimit. The intraband part of the susceptibility χ i in χ = χ i + χ e , (5)can be expanded in powers of ω Σ as ( − iω Σ ) χ i = η + ( − iω Σ ) σ + · · · , (6)or equivalently χ i = η z + σ z + · · · , (7)where z = − iω Σ . Clearly, χ i diverges at zero frequencysum. Since χ e is regular as ω Σ → χ itself divergesat zero frequency sum. As a side note, in metals , theFermi surface adds other divergent dc contributions tothe second-order BPVE, see Sec. VIII.Assuming a monocromatic optical field and usingMaxwell equation d P dt = J , (8)Eq. 6 implies η and σ are response functions of thenonlinear currents ddt J a (2) inj ≡ (cid:88) bc η abc (0 , ω, − ω ) E b ( ω ) E c ( − ω ) , (9) J a (2) sh ≡ (cid:88) bc σ abc (0 , ω, − ω ) E b ( ω ) E c ( − ω ) . (10) η and σ are the standard injection and shift cur-rent response functions derived from the susceptibil-ity approach. Importantly, they vanish for frequenciessmaller than the energy gap (they are ‘resonant’). Thedots in Eq. 6 are associated with the (nonresonant) rec-tification currents.
In the absence of momentum relaxation and saturationeffects the injection and shift currents grow with illumi-nation time as | J a (2) inj | ∝ η t, (11) | J a (2) sh | ∝ σ . (12)In this article we study how the injection and shift cur-rents are modified by the presence of a static field, i.e.,the fourth term in Eq. 1. We use the method of find-ing divergences of the free third-order electric polariza-tion susceptibility χ . Biased irradiated semiconductorsof this kind have been extensively studied numericallyusing the semiclassical Boltzmann equation. As shownbelow, this approach misses some important quantum ef-fects which are recovered in the susceptibility approach.In summary, the third order polarization P a (3) = (cid:88) bβcσdδ χ abcd ( − ω Σ , ω β , ω σ , ω δ ) E bβ E cσ E dδ e − iω Σ t , (13)oscillates with frequency ω Σ = ω β + ω σ + ω δ in the long-time limit. We show that the intraband part, χ i , of χ = χ i + χ e admits the Taylor expansion( − iω Σ ) χ i = ι + ( − iω Σ ) η + ( − iω Σ ) σ + · · · , (14)or alternatively the Laurent series χ i = ι z + η z + σ z + · · · , (15)where z = − iω Σ and ι , η , σ are (resonant) residues.Clearly, χ i diverges in the dc limit ( ω Σ = 0) and similarto η and σ , ι , η and σ represent response functionsof nonlinear currents d dt J a (3) jerk ≡ (cid:88) bcd ι abcd (0 , ω, − ω, E b ( ω ) E c ( − ω ) E d (16) ddt J a (3) inj ≡ (cid:88) bcd η abcd (0 , ω, − ω, E b ( ω ) E c ( − ω ) E d (17) J a (3) sh ≡ (cid:88) bcd σ abcd (0 , ω, − ω, E b ( ω ) E c ( − ω ) E d . (18)The difference is that a static field (zero frequency) istaken into account in addition to a monochromatic op-tical field. In the absence of momentum relaxation andsaturation effects the currents vary as t , t, t with illumi-nation time and we dub them jerk , third-order injectioncurrent and third-order shift current, respectively. Thedots in Eq. (15) represent regular terms associated withrectification currents.Since χ e is regular in the dc limit, one can write thesame expansion as in Eq. (15) for both χ i and χ . Sim-ilarly, the third order conductivity which is defined by J a (3) ≡ (cid:88) bβcσdδ σ abcd (3) ( − ω Σ , ω β , ω σ , ω δ ) E bβ E cσ E dδ e − iω Σ t , (19)admits the expansion σ (3) = ι z + η z + σ + · · · . (20)The subsequent evolution of charge distribution in thesample involves not only the above generation processesbut also the macroscopic current dynamics in a samplefor which momentum and energy relaxation is crucial.In the presence of dissipation, the dc divergences will becut off by a momentum relaxation time scale, just as thedc divergence of metals in the Drude model is cut off bya momentum relaxation time. In the BPVE, we expecttwo main relaxation time scales. One is the relaxationtime scale of the diagonal elements of the density ma-trix, τ , which ι , η , and η depend on. This could be ofthe order of 100 fs or longer in clean semiconductors. The second is the relaxation time scale of the off-diagonalelements of the density matrix, τ , which σ and σ de-pend on. Typically, τ < τ , but a recent experimentfound τ to be as large as 250 fs. For weakly disorderedsemiconductors, the photoconductivity ( ω Σ = 0) in Eq. 1becomes σ (3) ph ∼ τ ι + τ η + σ . (21)We can generalize the above results to any power in theelectric fields. In general, with each additional power inthe electric field, χ ni has an additional frequency factorin the denominator. This means that the dc singularitiesof χ ni are, at most, of the order n . We can show thatthe n th order z = 0 singularities of χ n ( n ≥ t n in the absence of mo-mentum relaxation and saturation effects. This occurswhen all but two of the external frequencies are zero.In addition, there is a hierarchy of higher order shift,injection,..., currents which are represented by z = 0 sin-gularities of order 1 , , , ..n of χ n . Formally χ n (or σ ( n ) )can be expanded as χ n = ∞ (cid:88) l = − n a l z l , (22)where a l = 0 for frequencies less than the gap and hencethe residues are a l = 12 πi (cid:73) | z | = ρ χ n dzz l +1 . (23)The poles of χ n may be of lower order than n when theoptical field is not monocromatic; see, for example, thefourth row in Table I where the field’s frequencies are ω and 2 ω .Importantly, we give simple physical arguments to ex-plain the microscopic processes involved in ι , η , σ and σ and provide explicit expressions in terms of ma-terial parameters amenable for first principles computa-tions. To have a sense of the magnitude of these cur-rents, we calculate them in single-layer GeS using a two-dimensional (2D) tight-binding model.The article is organized as follows. In Sec. II we de-scribe the conventions used in this paper. In Sec. III, IV,and V we introduce the Hamiltonian, polarization, andcurrent operators. In Sec. V A we revisit the calculationof the intraband current following Sipe and Shkrebtii. In Sec. VII we rederive the expressions for the injectionand shift current responses giving simple physical inter-pretations based on semiclassical wave packet dynamicsin electric fields. In Sec. VIII we include Fermi sur-face contributions to the second-order BPVE. We thenstudy the physical divergences of χ at zero frequency inSec. IX, X, XI, and XII. The jerk current has been pre-sented previously and is included here only for complete-ness . BPVEs arising from singularities of χ n ( n > II. NOTATION
To keep the notation under control we often omit theindependent variables such as time, real space position,or crystal momentum, specially in expressions which arediagonal in these variables. We use the standard notation for the n th electric po-larization susceptibility, χ abc...n ( − ω Σ , ω β , ω σ , ... ), where ω β , ω σ ,... label external frequency components, abc, ... label Cartesian components, and ω Σ = ω β + ω σ + ... thefrequency sum. We often write χ abc...n or simply χ n forbrevity absorbing a free permittivity factor (cid:15) into thesusceptibility.We adopt a semicolon and subscript, ‘ ; a ’, to mean a co-variant derivative with respect to crystal momenta withCartesian component a = x, y, z . Unless otherwise spec-ified, we contract spinor indices, e.g., nα → n in all ex-pressions. A hat on a Hamiltonian, polarization, andcurrent indicates an operator and a lack of a hat meansa quantum mechanical average. We do not use hats onthe creation or annihilation operators or on the positionoperator. A bold font indicates a vector or spinor.To distinguish the injection current derived from η from that of η we often call the former third-order in-jection current and the latter second-order injection cur-rent. Similarly, third-order shift current refers to currentderived from σ . We hope the missing details will becomeclear from the context. III. HAMILTONIAN
We start from a Hamiltonianˆ H = (cid:90) d r ψψψ † (cid:18) ˆ p m + V ( r ) + µ B e · (ˆ p × σσσ ) (cid:19) ψψψ, (24)describing Bloch electrons with spin-orbit (SO) coupling,where V ( r ) is the periodic potential of the ions, ˆ p = − i (cid:126) ∇∇∇ r is the momentum operator, e ( r ) = −∇∇∇ r V ( r ) isthe SO field from the nucleus, and µ B = e (cid:126) / mc is theBohr magneton. Electron-electron correlations in mean-field theory can be easily included by renormalizing theparameters of the noninteracting theory in Eq. 24. Mo-mentum relaxation is incorporated phenomenologicallyat the end of the calculation. The electron charge is e = −| e | . We define the real space spinor field as ψψψ = (cid:18) ψ ↑ ψ ↓ (cid:19) . (25)A classical homogeneous electric field is coupled to theHamiltonian by minimal substitution, ˆ p → ˆ p − e A . Afterthe gauge transformation˜ ψ α = ψ α e − ie A · r / (cid:126) , (26)( α is the spinor component) the Hamiltonian for thetransformed fields becomesˆ H ( t ) = ˆ H + ˆ H D ( t ) . (27)In what follows we omit the tilde above the transformedfields. ˆ H is given by Eq. (24), and the perturbation hasthe dipole formˆ H D = − e (cid:90) d r ψψψ † r · E ψψψ. (28)The electric field is given by E = − ∂ A /∂t . The eigen-functions of H can be chosen to be Bloch wavefunc-tions ψψψ ( β ) n ( kr ) = u ( β ) n ( kr ) e − i kr , where u ( β ) n ( k , r + R ) = u ( β ) n ( k , r ) has the period of a lattice vector R . k is thecrystal momentum and β = 1 , ψ α ( r ) = (cid:88) nβ k ψ ( β ) nα ( kr ) a nβ ( k ) , (29)where a † nβ ( k ) creates a particle in a Bloch state and obeysanticommutation rules { a † nα ( k ) , a mβ ( k (cid:48) ) } = δ nm δ αβ δ kk (cid:48) (= δ nm (2 π ) δ ( k − k (cid:48) ) /V in the thermodynamic limit).In this basis, H is diagonalˆ H = (cid:88) nβ k (cid:126) ω nβ a † nβ a nβ , (30)and (cid:126) ω nβ ( k ) is the energy of band n and spinor β . Thesum over crystal momenta is confined to the BrillouinZone (BZ). In the thermodynamic limit in d -dimensionsthe sum becomes (cid:80) k → V (cid:82) d d k/ (2 π ) d , where V is thevolume of the crystal. In what follows we chose the peri-odic gauge by which Bloch wavefunctions are periodic inreciprocal lattice vectors, ψψψ ( β ) n ( k + G , r ) = ψψψ ( β ) n ( k , r ). IV. POLARIZATION OPERATOR
The many-body polarization operator is well definedin finite systems. It is given byˆ P = 1 V (cid:90) d r ψψψ † e r ψψψ, (31)where e r /V is the one-body polarization operator. FromEq. 28, the dipole Hamiltonian becomes simplyˆ H D = − V ˆ P · E . (32)In periodic systems, H D is given in terms of Bloch oper-ators as ˆ P = eV (cid:88) nm kk (cid:48) (cid:104) n k | r | m k (cid:48) (cid:105) a † n ( k ) a m ( k (cid:48) ) . (33)Because the position operator is unbounded and theBloch wavefunctions extend to infinity, the matrix ele-ments (restoring spinor indices) (cid:104) n k | r | m k (cid:48) (cid:105) →(cid:104) nα k | r | mβ k (cid:48) (cid:105) = (cid:90) d r ψψψ ( α ) † n ( kr ) r ψψψ ( β ) m ( k (cid:48) r ) , (34) are singular. Fortunately, this singularity does not prop-agate to observables such as the spontaneous polariza-tion if we separate the singularity by the well-knownidentity (cid:104) n k | r | m k (cid:48) (cid:105) = δ nm [ δ ( k − k (cid:48) ) ξ nn + i ∇∇∇ k δ ( k − k (cid:48) )]+(1 − δ nm ) δ ( k − k (cid:48) ) ξ nm . (35)Here ξ nm are the Berry connections ξ nm → ξ nαmβ = (cid:90) d r u ( α ) † n i ∇∇∇ k u ( β ) m . (36)The polarization operator can then be separated into in-terband component proportional to (1 − δ nm ), and in-traband component proportional to δ nm . To tighten thenotation let us define the dipole matrix elements as r nm ≡ ξξξ nm n (cid:54) = m ≡ . (37)The polarization is then ˆ P = ˆ P e + ˆ P i , (38)where ˆ P e = eV (cid:88) nm k r nm a † n a m , (39)ˆ P bi = ieV (cid:88) n k a † n a n ; b , (40)and b = x, y, z . The intraband polarization depends onthe covariant derivative of a n a n ; b ≡ (cid:0) ∂∂k b − iξ bnn (cid:1) a n , (41)which transforms as a scalar, a n ; b → a n ; b e iφ βn , under localgauge transformations ψ ( β ) n → ψ ( β ) n e iφ ( β ) n . This should becontrasted with the transformation of ∂a n /∂k b which ac-quires a gauge-dependent contribution and hence it can-not represent a physical observable.From Eq. 38, the susceptibility also naturally separatesinto intraband and interband contributions as χ = χ i + χ e . (42) V. CURRENT OPERATOR
The current density is given byˆ J = eV (cid:90) d r ψψψ † ˆ v ψψψ, (43)where ˆ v = [ r , ˆ H ] /i (cid:126) = ˆ p /m + µ B σσσ × e is the electron’svelocity. In the presence of light, the momentum changesto ˆ p → ˆ p − e A , but after the gauge transformation (26),the current has the same expression. In terms of Blochoperators it becomesˆ J = eV (cid:88) nm k v nm a † n a m , (44)where v nm ≡ (cid:104) n k | ˆ v | m k (cid:105) . The current satisfies chargeconservation and Maxwell’s equation ∇ · ˆ j + ∂ ˆ ρ∂t = 0 (45) d ˆ P dt = ˆ J , (46)where ˆ ρ = eψψψ † ψψψ is the local charge density, ˆ j =( e/ ψψψ † ˆ v ψψψ +( e/ v ψψψ ) † ψψψ is the local charge current, andˆ P is the polarization given by Eq. 38. Local particle con-servation follows from the equation of motion (EOM) ofˆ ρ in the standard way. Maxwell’s equation is establishedas follows. From Eqs. 27 and 38 and i (cid:126) d ˆ P /dt = [ ˆ P , ˆ H ],we obtain i d ˆ P a dt = eV (cid:88) nm k (cid:0) iω n ; a δ nm + ω mn r anm (cid:1) a † n a m (47)where ω nm ≡ ω n − ω m . We define the covariant derivativeof the matrix element O nm ≡ (cid:104) n k | O | m k (cid:105) between Blochstates n, m at a single crystal momentum by O nm ; b ≡ (cid:20) ∂∂k b − i ( ξ bnn − ξ bmm ) (cid:21) O nm , (48)which can be shown to transform as a tensor under gaugetransformations. Since the energy bands are the diag-onal matrix elements of the Hamiltonian, their covari-ant derivative reduces to the standard derivative ω n ; a = ∂ω n /∂k a = v an = p an /m + µ B ( σσσ × e ) ann . On the right handside of Eq. 47, we recognize the diagonal and off-diagonalmatrix elements of the velocity. The off-diagonal matrixelements are obtained by taking Bloch matrix elementson both sides of ˆ v = [ r , ˆ H ] /i (cid:126) . Comparing with Eq. 44,the Maxwell’s equation is established in the basis of Blochoperators.The intraband polarization operator defines the intra-band current operator which, as shown below, connectsthe semiclassical wave packet dynamics and the BPVEs. A. Intraband current
We define the intraband current operator as the timederivative of the intraband polarization operator ˆ J i = d ˆ P i /dt . Similarly, the interband current is ˆ J e = d ˆ P e /dt .The total current is the sum of the twoˆ J = ˆ J i + ˆ J e . (49) Let us first calculate ˆ J e from i (cid:126) d ˆ P ae dt = [ ˆ P ae , ˆ H ] − V (cid:88) b [ ˆ P ae , ˆ P bi + ˆ P be ] E b . (50)The first term has been computed in Eq. (47). The sec-ond term is[ ˆ P ae , ˆ P bi + ˆ P be ] = − ie V (cid:88) nm k (cid:0) r anm ; b + i (cid:88) p [ r anp r bpm − r bnp r apm ] (cid:1) a † n a m . (51)To make progress we now invoke a sum rule first discussedby Sipe and coworkers. It derives from taking matrixelements of [ r a , r b ] = 0 , (52)and carefully separating the interband and intrabandparts of the position operator shown in Eq. 35. It is easyto show that such procedure works for spinor matrix el-ements too. Two cases are of interest follow. Takingdiagonal matrix elements ( n = m ) of Eq. 52 givesΩ ban ≡ ∂ξ ann ∂k b − ∂ξ bnn ∂k a = − i (cid:88) l [ r anl r bln − r bnl r aln ] , (53)and off-diagonal elements ( m (cid:54) = n ) gives r anm ; b − r bnm ; a = − i (cid:88) l [ r anl r blm − r bnl r alm ] . (54)It is customary, in analogy with electrodynamics, to de-fine a gauge field tensor Ω abn derived from the Berry vectorpotential of band n . The Berry curvature ΩΩΩ n = ∇∇∇ × ξξξ nn is related to the gauge field by Ω abn = (cid:80) e (cid:15) abe Ω en . Wenow separate the diagonal from the nondiagonal matrixelements in Eq. 51 and use Eqs. (53,54) to obtain − V (cid:88) b [ ˆ P ae , ˆ P bi + ˆ P be ] E b = ie V (cid:88) n k ( E × ΩΩΩ n ) a a † n a n + ie V (cid:88) nm k b E b r bnm ; a a † n a m . (55)Subtracting ˆ J e (Eq. 50) from ˆ J (Eq. 47) we obtain ˆ J i ˆ J ai = eV (cid:88) nm k (cid:20) ω n ; a δ nm − e (cid:126) ( E × ΩΩΩ n ) a δ nm − e (cid:126) E · r nm ; a (cid:21) a † n a m . (56)This is an important result. The first term is the stan-dard group velocity (renormalized by the SOC) of anelectron wave packet in band n , ω n ; a = v an . As shownbelow, this term gives rise to the injection current con-tribution to the BPVE. The second term depends on theBerry curvature ΩΩΩ n and is often called ‘anomalous’ ve-locity. It gives rise to many topological effects in con-densed matter physics. For example, it gives rise tothe (intrinsic) anomalous Hall conductivity in metallicferromagnets, and, as shown in Sec. VIII, to the (in-trinsic) nonlinear Hall effect in nonmagnetic metals. In insulators, this term contributes to third order in theelectric field but not to second order.The third term resembles a small dipole created bythe external electric field. Just as the standard momen-tum derivative of Bloch energies leads to the usual groupvelocity, the (covariant) derivative of the dipole energy U nm = e E · r nm , can be thought of as a group velocity v adip,nm = − e (cid:126) E · r nm ; a (57)associated with a pair of wave packets in distinct bands.The first two integrands in Eq.(56) are gauge invariantand are usually interpreted as velocity contributions ofelectron wave packets. The dipole velocity, on the otherhand, is not gauge invariant and hence is not a physicalvelocity. However, the product of the dipole velocity andthe density matrix is gauge invariant and in this contextthe dipole velocity can be given the interpretation of thevelocity of pairs of wave packets. As shown below, thedipole velocity gives rise to the shift current contribu-tion to the BPVE. Indeed, the intraband current unifiesthe well-known semiclassical dynamics of wave packets inelectric fields with the BPVEs.What is the physical interpretation of the interbandcurrent? The fact that the interband polarization is regu-lar in the dc limit ( ω Σ →
0) implies the interband currentvanishes in this limit. This suggests that the interbandcurrent captures electron oscillations about their equilib-rium positions but not their uniform acceleration.Up to this point, the above formalism is valid for met-als and insulators. Except for Sec. VIII, we will fo-cus on the short time response of insulators, discardingFermi surface contributions and momentum relaxation.By ‘short time’ we mean shorter than momentum relax-ation characteristic time ( ∼
100 fs) but longer than theperiod of light ( ∼ VI. PERTURBATION THEORY
Let us define the single-particle density matrix ρ mn ≡ (cid:104) a † n a m (cid:105) , (58)where the a n operators are in the Heisenberg represen-tation. The quantum average is over the ground statedefined with all the valence bands filled and all conduc-tion bands empty. Being noninteracting, the system is completely characterized by the single-particle densitymatrix. The amplitude of the electric field is E b = (cid:88) β E bβ e − i ( ω β + i(cid:15) ) t , (59)where β = 1 , , ... labels the frequency components of thefield. The dipole Hamiltonian is treated as a perturba-tion with the electric field being turned on slowly in theinfinite past so that all the transients effects have van-ished. As usual, this is accomplished by taking the limit (cid:15) → ∂ρ mn ∂t + iω mn ρ mn = ei (cid:126) (cid:88) lb E b ( ρ ml r bln − r bml ρ ln ) − e (cid:126) (cid:88) b E b ρ mn ; b . (60)The first term on the right comes from interband pro-cesses as can be recognized by the presence of r nm . Thesecond term comes from intraband processes which in-volve the covariant derivative of the density matrix ρ mn ; b ≡ (cid:20) ∂∂k b − i ( ξ bmm − ξ bnn ) (cid:21) ρ mn , (61)Only when the intraband and interband motion is consid-ered on an equal footing, the EOM reduces to the Boltz-mann equation (in the one-band limit) with no collisionintegral. A. Zeroth order If E = 0 the solution of Eq.(60) is simply ρ (0) mn = δ nm f n ,where f n ≡ f ( (cid:15) n ( k )) = 0 , n at zero temperature. B. First order
Substituting the zeroth-order solution into the right-hand side of Eq.(60) and solving for ρ (1) mn we obtain ρ (1) mn = (cid:88) bβ ¯ ρ (1) bβmn E bβ e − iω β t (62)= e (cid:126) (cid:88) bβ r bmn f nm ω mn − ω β E bβ e − iω β t . (63)where we defined f nm ≡ f n − f m . Note that to first orderonly interband processes are allowed in insulators. C. Second order
To second order we have ρ (2) mn = (cid:88) bβ (cid:88) cσ ¯ ρ (2) bβcσmn E bβ E cσ e − iω Σ t , (64)where¯ ρ (2) bβcσmn = ie (cid:126) ( ω mn − ω Σ ) (cid:20) ¯ ρ (1) bβmn ; c + i (cid:88) l (cid:0) ¯ ρ (1) bβml r cnl − r cml ¯ ρ (1) bβln (cid:1)(cid:21) , (65)and ω Σ = ω β + ω σ . The covariant derivative of a quotientin ¯ ρ (1) bβmn ; c is simply (cid:18) r amn f nm ω mn − ω α (cid:19) ; b = r amn ; b f nm ω mn − ω α − r amn f nm ω mn ; b ( ω mn − ω α ) (66) D. nth-order
In the long-time limit, we expect harmonic solutionsof the form ρ ( n ) mn = (cid:88) a α ,... ¯ ρ ( n ) a α ,..mn E a α · · · E a n α n e − iω ( n )Σ t , (67)where ω ( n )Σ = ω α + · · · + ω α n . Substituting into Eq.(60)and iterating we obtain an equation for ¯ ρ ( n +1) mn in termsof ¯ ρ ( n ) mn . Omitting the supercripts a α , ... for clarity weobtain¯ ρ ( n +1) mn = ie (cid:126) ( ω mn − ω ( n +1)Σ ) (cid:20) i (cid:88) l (¯ ρ ( n ) ml r a n +1 ln − r a n +1 ml ¯ ρ ( n ) ln )+ ¯ ρ ( n ) mn ; a n +1 (cid:21) . (68)Note that at every order in perturbation theory thereare interband (first term) and intraband (second term)contributions. In general, the n th-order ρ ( n ) ( n ≥
1) has2 n − intraband and 2 n − interband contributions. VII. PHYSICAL DIVERGENCES OF χ The susceptibility and conductivity response tensorsto second order are defined by P a (2) = (cid:88) bβcσ χ abc ( − ω Σ , ω β , ω σ ) E bβ E cβ e − iω Σ t , (69) J a (2) = (cid:88) bβcσ σ abc (2) ( − ω Σ , ω β , ω σ ) E bβ E cβ e − iω Σ t , (70) where ω Σ = ω β + ω σ . They are related by dP a (2) /dt = J a (2) . χ can be split into interband and intraband com-ponents, χ = χ e + χ i , using Eqs.(39),(56),(63), and(64). The result is χ abc e C = i (cid:88) nm k r anm f nm ω mn − ω Σ (cid:18) r bmn ω mn − ω β (cid:19) ; c − (cid:88) nlm k r anm ω mn − ω Σ (cid:18) r bml r cln f lm ω ml − ω β − r cml r bln f nl ω ln − ω β (cid:19) , (71) χ abc i C = iω (cid:88) nm k ω nm ; a r bnm r cmn f mn ω nm − ω β + 1 iω Σ (cid:88) nm k r cnm ; a r bmn f nm ω mn − ω β , (72)where we defined C = e / (cid:126) V . These expressions needto be symmetrized with respect to exchange of indices bβ ↔ cσ . We note that χ i is easier to calculate from J (2) i rather than directly from P (2) i .The Taylor expansion of χ i in Eq. 6 means that χ i diverges as ω Σ → η and shift σ response tensors can be obtained from this expansion,see Appendix B. Here we derive these tensors from aslightly different perspective that exposes the analyticproperties of χ i . Let us assume χ i admits a Laurentseries χ i = η z + σ z + · · · (73)where z = − iω Σ . Then η is given by η = 12 πi (cid:73) | z | = ρ dz zχ i , (74)where ρ is the radius of convergence. All the frequen-cies are parametrized in terms of ω Σ = iz . One suchparametrization is ω β = ω + n β ω Σ (75) ω σ = − ω + n σ ω Σ , (76)where n β + n σ = 1. The manifold where ω Σ = 0 is a lineof singular points ( ω β , ω σ ) = ( ω, − ω ), parametrized by asingle frequency ω >
0. Symmetrizing χ i with respect toexchange of indices bβ ↔ cσ and using Eq. 74 we obtain η abc (0 , ω, − ω ) as η abc = πe (cid:126) V (cid:88) nm k f mn ω nm ; a r bnm r cmn δ ( ω nm − ω ) , (77)or equivalently η abc = πe (cid:126) V (cid:88) nm k f mn ω nm ; a ( r bnm r cmn − r cnm r bmn ) × δ ( ω nm − ω ) , (78)which is independent of the parameters n β , n σ . In calcu-lating η we take the limit ρ → (cid:15) → ω Σ = 0in the infinite past. Similarly, σ abc (0 , ω, − ω ) is given by σ = 12 πi (cid:73) | z | = ρ dz χ i . (79)An explicit integration gives σ abc = iπe (cid:126) V (cid:88) nm k f mn ( r cnm ; a r bmn − r cnm r bmn ; a ) δ ( ω nm − ω ) . (80)In calculating σ we took n β = n σ = 1 / n β − n σ . Thisterm does not arise in the standard method becausethere the prescription is to Taylor expand only the realparts. Taking n β = n σ means we are approaching theline of singularities at right angle.Eqs.(77) and (80) are the well-known injection andshift current tensors. η is pure imaginary and antisym-metric in the b, c indices and hence vanishes for linearpolarization. σ , on the other hand, is real, symmetric in b, c indices and hence vanishes for circular polarization.The corresponding injection and shift currents are givenby J a (2) sh ≡ (cid:88) bβcσ σ abc ( − ω Σ , ω β , ω σ ) E bβ E cσ e − iω Σ t , (81) ddt J a (2) inj ≡ (cid:88) bβcσ η abc ( − ω Σ , ω β , ω σ ) E bβ E cσ e − iω Σ t , (82)subject to ω Σ = 0. Assuming a monocromatic source E b = E b ( ω ) e − iωt + c.c. and performing the frequencysums keeping only dc terms ( ω Σ = 0), we obtain J a (2) sh = 2 (cid:88) bc σ abc (0 , ω, − ω ) E b ( ω ) E c ( − ω ) (83) ddt J a (2) inj = 2 (cid:88) bc η abc (0 , ω, − ω ) E b ( ω ) E c ( − ω ) , (84)where the factor of 2 is from the intrinsic permutationsymmetry of susceptibilities. Being quadratic in thefields the injection and shift current vanish for centrosym-metric systems. The above expressions indicate the in-jection and shift currents vary as | J (2) inj ( t ) | ∼ η t (85) | J (2) sh ( t ) | ∼ σ (86)with illumination time in the absence of momentum re-laxation and saturation effects. A. Physical interpretation of injection and shiftcurrent
In this section we show that the injection and shift cur-rents can be understood from simple semiclassical wavepacket dynamics in electric fields.
1. Injection current
The microscopic origin of the injection current fromlight-matter interactions is well known. It arises from theasymmetry in the carrier injection rate at time-reversedmomenta in the BZ . To see this, let us consider anelectron wave packet with velocity v an . From the firstterm in Eq. 56 the current is J a = eV (cid:88) n k f n v an , (87)where f n ≡ ρ (0) nn . The effect of an optical field is to inject carriers into the current-carrying states in the conductionbands. Taking a time derivative of the occupations weobtain ddt J ainj = eV (cid:88) n k df n dt v an . (88)For low intensity, the Fermi’s Golden Rule gives the one-photon absorption rate df v dt = − πe (cid:126) (cid:88) c | E ( ω ) · r cv | δ ( ω cv − ω ) ,df c dt = 2 πe (cid:126) (cid:88) v | E ( ω ) · r cv | δ ( ω cv − ω ) , (89)where c, v labels a conduction or a valence band, respec-tively. For complex fields, e.g, circularly polarized orelliptically polarized, the carrier injection rate at time-reversed points ± k in the BZ is not the same ddt f c ( − k ) (cid:54) = ddt f c ( k ) , (90)leading to a polar distribution of Bloch velocity states.This is the microscopic origin of the injection currentand, as we show below, of many higher-order injectioncurrents. Substituting into Eq.(88) we obtain ddt J a (2) inj = 2 πe (cid:126) V (cid:88) b (cid:48) c (cid:48) (cid:88) cv k ω cv ; a r b (cid:48) vc r c (cid:48) cv δ ( ω cv − ω ) × E b (cid:48) ( ω ) E c (cid:48) ( − ω ) , (91)or ddt J a (2) inj = 2 πe (cid:126) V (cid:88) bc (cid:88) nm k f mn ω nm ; a r bnm r cmn δ ( ω nm − ω ) × E b ( ω ) E c ( − ω ) , (92)0 � m e spaceE x (t) ω E x (t) ρ (t)momentum ρ (t) out of phase in-phaseconduction bandvalence band currentx FIG. 1. Intuitive picture of microscopic generation of shiftcurrent. The wiggle lines represent interband coherence os-cillations between the valence and conduction band centersof charge (circles) which are spatially separated. The quan-tum interference between population oscillations ρ nm ( t ) anddipole velocity oscillations E ( t ) · r mn ; a gives rise to a shiftcurrent. which is the standard injection current shown in Eq.(84).
2. Shift current
Injection current is proportional to the momentum re-laxation time and hence explicitly breaks time-reversalsymmetry. In the scenario where a shift current origi-nates from light-matter interactions, the shift currentdoes not require the presence of momentum relaxationto break time reversal symmetry. How is time-reversalsymmetry broken in shift current processes? It is brokenat the time of photon absorption which is an irreversibleprocess.Materials that exhibit shift current have valence andconduction band centers spatially separated within theunit cell and hence charge is shifted upon photon absorp-tion. This process depends only on the off-diagonal ele-ments of the density matrix and hence it requires quan-tum coherence as has been extensively documented. Herewe propose that shift current arises from the quantuminterference of two distinct microscopic processes involv-ing wave packet oscillations in the presence of an electricfield. To see this consider the dipole current in Eq.(56)to second order J a (2) dip = − e (cid:126) V (cid:88) nm k E ( t ) · r nm ; a ρ (1) mn ( t ) . (93)The current is the sum of dipole velocities of each pair ofwave packets in bands n, m weighted by the probability ρ (1) mn of being occupied. From Eq. 63 we have J a (2) dip = − e (cid:126) V (cid:88) bβcσ (cid:88) nm k r bnm ; a r cmn f nm ω mn − ω σ E bβ E cσ e − iω Σ t , (94)where ω Σ = ω β + ω σ . Symmetrizing with respect toexchange of indices bβ ↔ cσ , assuming a monocromaticfield E b = E b ( ω ) e − iωt + c.c. , and keeping only the dc resonant terms we obtain J a (2) sh = iπe (cid:126) V (cid:88) bc (cid:88) nm k f mn ( r bnm ; a r cmn + r cnm ; a r bmn ) × δ ( ω nm − ω ) E b ( ω ) E c ( − ω ) , (95)which is the standard expression for the shift current inEq. 83. This calculation suggests that the constructive quantum mechanical interference of interband coherenceoscillations and dipole velocity oscillations is the micro-scopic origin of the shift current, see Fig. 1. We note thatelectron oscillations between centers of charge, alone, donot lead to a dc current. However, the directionality ofthe electron oscillations combined with an isotropic re-laxation (due to, e.g., randomized collisions) could, inprinciple, also lead to a dc current. In this scenario mo-mentum relaxation would play a significant role in theorigin of the current.Before showing how the injection and shift currentsare modified by the presence of a static electric field,we discuss Fermi surface contributions to the BPVE tosecond order. VIII. THE BPVE IN METALS
The Fermi surface of metals gives rise to two additionalcontributions to the second-order BPVE. The first is thenonlinear Hall effect (NLHE) discussed by Sodemann andFu and the second is a metallic jerk current discussedrecently by Matsyshyn and Sodemann . Here we showthat these photocurrents can be obtained from Eq. 56and simple physical assumptions.To begin note that ρ (1) nm has a Fermi surface contribu-tion. Substituting ρ (0) nm = f n δ nm into the right-hand sideof Eq. 60 two terms are obtained. The first is an inter-band contribution given by Eq. 63. The second is theintraband contribution ρ (1) nm,i = − δ nm ie (cid:126) (cid:88) bβ ω β ∂f n ∂k b E bβ e − iω β t , (96)which depends explicitly on the presence of a Fermi sur-face via ∂f n /∂k b . The dc divergence is cut off by themomentum relaxation time scale τ as1 − iω β → τ − iω β . (97)From the second term in Eq. 56 we have J a (2) nlhe = − e (cid:126) V (cid:88) bβ (cid:88) n k e (cid:15) abe E bβ e − iω β ΩΩΩ en ρ (1) nn,i , (98)after symmetrizing with respect to exchanges of field in-dices, and performing the frequency sums the dc currentis1 J a (2) dc,nlhe = (cid:88) bc σ abc (2) nlhe (0 , ω, − ω ) E b ( ω ) E c ( − ω ) + c.c. (99)where σ abc (2) nlhe (0 , ω, − ω ) ≡ e (cid:126) V τ iωτ (cid:88) n k e (cid:15) abe ΩΩΩ en ∂f n ∂k c , (100)is the known response tensor for the dc nonlinear Halleffect. From a semiclassical point of view the dc NLHEarises from the quantum interference of (intraband) oscil-lations of excitations across the Fermi surface and oscil-lations of the anomalous velocity of wave packets. A keydifference with injection current is that NLHE responsetensor is antisymmetric in the first two indices, ratherthan the last two.Similarly, inspection of Eq. 56 shows that another dccurrent to second-order is possible by taking two timederivatives of the velocity in the first term and the equi-librium density matrix d J a dt = eV (cid:88) n k f n d v a dt . (101)The two derivatives of the velocity can be computedin powers of the electric field. The result is similar toEq. 126 but with the optical fields replacing the staticfields. Symmetrizing and performing the frequency sumswe obtain J a (2)2 ndjerk = (cid:88) bc ι abc (0 , ω, − ω ) E b ( ω ) E c ( − ω ) (102)where ι abc (0 , ω, − ω ) ≡ e (cid:126) V (cid:88) n k f n ω n ; abc . (103)From the semiclassical perspective, the second-order jerkcurrent arises from the constant acceleration of a wavepacket. The acceleration can be constant because of the(classical) interference of the oscillating electric field withitself. The response tensor is symmetric under exchangeof b, c indices and hence vanish for circular polarization.The second-order jerk current varies as t in the absenceof momentum dissipation and saturation effects.The responses at 2 ω can be calculated similarly. Theresults are J a (2)2 ω,nlhe = (cid:88) bc σ abc (2) nlhe ( − ω, ω, ω ) E b ( − ω ) E c ( − ω ) + c.c., (104) TABLE II. Second-order BPVEs in metals. intra=intraband,inter=interband, FS= Fermi surface. ω Σ = ω β + ω σ , frequencysum of two optical fields, ∗ in the absence of momentum re-laxation.dc Symbol Time intra. vs Sing. ∗ Mom.current dep. ∗ inter. vs relax.FSInjection η t intra ω Σ = 0 τ Shift σ const. inter ω Σ = 0 τ NLHE σ (2) nlhe const. intra,FS ω β = 0 τ Jerk ι t intra,FS τ where σ abc (2) nlhe ( − ω, ω, ω ) = σ abc (2) nlhe (0 , ω, − ω ) and J a (2)2 ndjerk = (cid:88) bc ι abc ( − ω, ω, ω ) E b ( ω ) E c ( − ω ) + c.c. (105)where ι abc ( − ω, ω, ω ) = (1 / ι abc (0 , ω, − ω ). Note thatboth, the NLHE and the second-order jerk current canproduce current transverse to the polarization of the opti-cal field. This should be taken into account in interpret-ing experiments in metals. A summary of the metallicBPVEs is given in Table II. IX. PHYSICAL DIVERGENCES OF χ The susceptibility and conductivity response tensorsto third order are defined by P a (3) = (cid:88) bβcσdδ χ abcd ( − ω Σ , ω β , ω σ , ω δ ) E bβ E cσ E dδ e − iω Σ t , (106) J a (3) = (cid:88) bβcσdδ σ abcd (3) ( − ω Σ , ω β , ω σ , ω δ ) E bβ E cσ E dδ e − iω Σ t , (107)where ω Σ = ω β + ω σ + ω δ . They are related by d P (3) /dt = J (3) . χ can be split into interband and intraband com-ponents χ = χ e + χ i . Expanding the intraband com-ponent in powers of ω Σ gives( − iω Σ ) χ i = ι + ( − iω Σ ) η + ( − iω Σ ) σ + · · · . (108)See Appendix C. Eq.(108) is equivalent to χ i = ι z + η z + σ z + · · · (109)where z ≡ − iω Σ . Since χ e is regular, Eq. 109 impliesthat the conductivity in the limit of no momentum re-laxation is σ (3) = ι z + η z + σ + z ( reg ) , (110)2where reg represents the remaining regular terms (as z → ι , η and σ define various currentcontributions as follows. The limitlim z → z σ (3) = ι , (111)or equivalentlylim ω Σ → d dt J a (3) ≡ d dt J a (3) jerk = (cid:88) bβcσdδ ι abcd (0 , ω β , ω σ , ω δ ) E bβ E cβ E dδ , (112)(subject to ω Σ = 0) defines the jerk current. Similarlythe limits lim z → z (cid:104) σ (3) − ι z (cid:105) = η , (113)lim z → (cid:104) σ (3) − ι z − η z (cid:105) = σ , (114)define higher order injection and shift currents (respec-tively) in the presence of a static electric field: ddt J a (3) inj ≡ (cid:88) bβcσdδ η abcd (0 , ω β , ω σ , ω δ ) E bβ E cβ E dδ , (115) J a (3) sh ≡ (cid:88) bβcσdδ σ abcd (0 , ω β , ω σ , ω δ ) E bβ E cβ E dδ , (116)subject to ω Σ = 0. We now analyze each of these currentsin detail. X. JERK CURRENTA. Hydrodynamic model
In an isotropic system the current is J aclas = env a , (117)where n is the carrier density. Taking two derivatives weobtain d dt J aclas = e d ndt v a + 2 e dndt dv a dt + en d v a dt . (118)If the rate of carrier injection dn/dt = g and acceleration eE a /m ∗ are constant in time then d dt J aclas = 2 e gE a m ∗ = constant, (119)leads to a current varying quadratically with illuminationtime. This effect has been extensively studied in the con-text of the THz generation in bias semiconductor anten-nas using semiclassical kinetic equations, see for exampleRef. 52. However, the static field modifies the carrier in-jection rate giving rise to novel contributions. We nowdiscuss this effect. B. Susceptibility divergence
We find ι from the limit lim ω Σ → ( − iω Σ ) χ i = ι .The details of the derivation are outlined in Appendix D. ι abcd (0 , ω, − ω,
0) is given by ι abcd = 2 πe (cid:126) V (cid:88) nm k f mn (cid:2) ω nm ; ad r bnm r cmn + ω nm ; a ( r bnm r cmn ) ; d (cid:3) δ ( ω nm − ω ) , (120)where ω nm ; ad = ∂ ω nm /∂k d ∂k a = ∂ ω n /∂k d ∂k a − ∂ ω m /∂k d ∂k a .Assuming time-reversal symmetry in the groundstate we can choose r nm ( − k ) = r mn ( k ) to showthat ι is real, symmetric in the b, c indices, andsatisfies [ ι abcd (0 , ω, − ω, ∗ = ι acbd (0 , ω, − ω,
0) = ι abcd (0 , − ω, ω, ι controlsthe current d dt J a (3) jerk = (cid:88) bβcγdδ ι abcd ( − ω Σ , ω β , ω γ , ω δ ) E bβ E cγ E dδ e − iω Σ t , (121)subject to ω Σ = 0. Performing the sum over frequencieswe obtain d dt J a (3) jerk = 6 (cid:88) bcd ι abcd (0 , ω, − ω, E b ( ω ) E c ( − ω ) E d , (122)where E d is a static external field. The factor of 6 = 3!is the number of pair-wise exchanges of field indices( bβ ) , ( cσ ) , ( dδ ). The jerk current vanishes for frequen-cies smaller than the energy band gap. Eq.(122) indicatesthat the jerk current grows quadratically with illumina-tion time | J (3) jerk ( t ) | ∼ ι t , (123)in the absence of momentum relaxation and saturation ef-fects. In analogy with second derivative of velocity whichis called ’jerk’ we dub it jerk current. This should becompared and contrasted with injection current whichgrows linearly with illumination time (Eq. 85) and shiftcurrent which is constant (Eq. 86). C. Materials
In general, the 81 components of ι are finite in bothcentrosymmetric and noncentrosymmetry crystal struc-tures. In practice, the symmetries of the 32 crystal classesgreatly reduce the number of independent components.For example, GaAs has ¯43 m point group, with 21 nonzerocomponents and 4 independent components. However, ι is symmetric under exchange of bc which reduces thenumber of independent component to 3. In 2D materi-als the number of components of ι is also small. For3example, single-layer GeS has mm ι has only six independent components.In general, linear, circular or unpolarized light will pro-duce jerk current along the direction of the static field.Current transverse to the static field may not be gener-ated with unpolarized or circular polarization. D. Physical interpretation of jerk current
The terms in Eq.(120) are hard to interpret physically.We now rederive the same result in a physically moretransparent way using a phenomenological model. Con-sider an electron wave packet in band n subject to a staticelectric field E d . The electron’s wavevector obeys (cid:126) d k dt = − e ∂ A ∂t , (124)where the vector potential A gives the static electricfield E d = − ∂A d /∂t . The Bloch velocity of the electron v n ( k − e A / (cid:126) ) can be expanded in powers of A . Its firstand second time derivatives are given by dv an dt = e (cid:126) (cid:88) d ω n ; ad E d , (125) d v an dt = e (cid:126) (cid:88) de ω n ; ade E d E e . (126)Now, taking two time derivatives of Eq.(87) d J a dt = eV (cid:88) n k (cid:18) d f n dt v an + 2 df n dt dv an dt + f n d v an dt (cid:19) , (127)and using Eq. 89 and Eq. 125 we have (to linear order in E d ) d dt J a (3) jerk =2 πe (cid:126) V (cid:88) bc (cid:48) d (cid:88) cv k ω cv ; ad r bvc r c (cid:48) cv δ ( ω cv − ω ) E b ( ω ) E c (cid:48) ( − ω ) E d + 2 πe (cid:126) V (cid:88) bc (cid:48) d (cid:88) cv k ω cv ; a ∂ ( r bvc r c (cid:48) cv ) ∂k d δ ( ω cv − ω ) × E b ( ω ) E c (cid:48) ( − ω ) E d . (128)Since ω > d f n /dt (cid:54) = 0 whichis missing in the standard semiclassical approach. - - E ne r g y ( e V ) GeS Γ X S Y Γ (a) xy t t t’ t’ t a Top viewSide view P XM Z d (b) y x xy Tight-binding V Ld Setup li gh t (c) xy L FIG. 2. (a) Band structure of single-layer GeS indicat-ing transitions near the band edge (red arrow). (b) crys-tal structure of single-layer GeS, (c) sample setup and two-dimensional, two-band tight binding model of single-layer GeSwhich reproduces the nonlinear optical response of near theband edge. The hopping parameters considered are indicated.See main text for more details.
E. Jerk Hall current
In an isotropic medium, charge carriers move paral-lel to the electric field. The jerk current, on the otherhand, can flow transverse to the static electric field ina rotationally symmetric medium. To see this, let usassume a sample biased in the x -direction and com-pute the current in the y -direction. An optical field E = ˆ x E x ( ω ) e − iωt + ˆ y E y ( ω ) e − iωt + c.c. , with E a ( ω ) = | E a ( ω ) | e − iφ a , is incident perpendicular to the samplesurface which defines the xy -plane. The current in the y − direction is4 d dt J y (3) jerk = ς yx jH E x , (129)where the effective jerk Hall (jH) conductivity is ς yx jH ≡ ι yxxx | E x ( ω ) | + 6 ι yyyx | E y ( ω ) | + 12 ι yyxx | E x ( ω ) || E y ( ω ) | cos( φ x − φ y ) . (130)In a simple relaxation time approximation the jerk con-ductivity (see Eq. 110) is cut off by a relaxation time τ as ι ( − iω Σ ) → ι ( τ − iω Σ ) (131)Hence, the jerk current is J y (3) jerk ∼ τ (1 − iω Σ τ ) ς yx jH E x , (132)where τ is the relaxation time of the diagonal elementsof the density matrix. In the dc limit the current is pro-portional to the square of the momentum relaxation. Forfrequencies larger than the Drude peak but smaller thaninterband transitions the current is independent of thescattering time and it is a measure of the geometry ofthe Bloch wavefunctions.The dependence on light’s polarization as cos( φ x − φ y )and the square of the momentum relaxation are uniquecharacteristics of the jerk current which can be used todistinguish it from η and σ .The symmetries of the crystal can also constrain thecontributions to the jerk current, e.g., if the crystal hasmirror symmetry y → − y the first and second termsin Eq. 130 vanish. In addition, for circular polarization φ x − φ y = ± π/ F. Example: Jerk current in single-layer GeS
To get a sense of the magnitude of the jerk current inreal materials we now calculate it for single-layer GeS.Single-layer GeS is of great interest for its predicted in-plane spontaneous ferroelectric polarization, suitable en-ergy band gap in the visible spectrum ( ∼ We consider a 2D, two-band tight-binding model ofsingle-layer GeS shown in Fig. 2c. The details of themodel are presented in Appendix G. The model has beenshown to reproduce the ab-initio shift and injection cur-rent of single-layer GeS near the band edge, specif-ically in the energy range 1.9-2.14 eV. Since the modelis 2D, we divide the model’s 2D current by the thicknessof the GeS layer ( d ∼ . ι ab c d ( X A m / V s ) xxxxxxyy (x10) xyyxyxxyyyxx (x10) yyyy Photon energy (eV)
FIG. 3. Jerk current response tensor of single-layer GeS nearthe band edge. The two-band model used is shown in Fig. 2c.The tensor vanishes for photon energies lower than the energyband gap ( ∼ . ). The strongest component is alongthe polar axis xxxx . The components yyxx, xxyy , describe aHall-like response and are an order of magnitude smaller. Foradded clarity, these components are multiplied by 10. Because of the mirror symmetry y → − y of the crystal,only six tensor component are independent. As seen inFig. 3, the strongest is along the polar (chosen along x -axis) of magnitude ∼ Am/V s . The current trans-verse to the static electric field, described by the compo-nent ι yyxx (see Eq. 129), is an order of magnitude smaller.The sample is rectangular of dimensions L × L andthickness d = 2 .
56 ˚A and is biased by an external batteryof voltage V , as seen Fig. 2c. Let us assume the opticalfield is incident perpendicularly to the plane of single-layer GeS as E ( t ) = ˆ x E x ( ω ) e − iωt + ˆ y E y ( ω ) e − iωt + c.c. (133) E = ˆ x E x . (134)where E x ( ω ) = E ( ω ) cos θe − iφ x , E y ( ω ) = E ( ω ) sin θe − iφ y , θ is the angle with the x -axis.The longitudinal and transverse currents are I x (3) jerk = 6 Aτ ( ι xxxx | E x ( ω ) | + ι xyyx | E y ( ω ) | ) E x , (135) I y (3) jerk = 12 Aτ ι yyxx | E x ( ω ) || E y ( ω ) | cos( φ x − φ y ) E x , (136)where A = Ld is the transverse area of the sample. Notethat the current along the polar ( x )-axis is independentof the polarization of light. Hence, the polar componentof the current will not vanish even for unpolarized light.The transverse component of the current, on the otherhand, vanishes for circularly polarized (and unpolarized)light and is maximum for linearly polarized light.We choose the optical field to be linearly polarized( φ x = φ y ) at an angle θ with the polar axis as shown inthe inset to Fig. 4a. The figure shows the jerk current in-duced as a function of θ . We assumed semiconductor pa-rameters typically found in the laboratory: L = 100 µ m,5 Photon energy (eV) I ( n A ) I ( p A ) θ=0 V E E(t) θ xy (a)(b) I x I y xy FIG. 4. Jerk current in single-layer GeS with linear polar-ization at various angles θ with respect to the polar axis.(a) current parallel to the polar axis and (b) current perpen-dicular to the polar axis. The transverse current is largestat θ = 45 whereas the parallel current is largest when thelight’s polarization is along the polar axis. The inset showsthe top view of the sample. V = 1V, E x = V /L = 10 V/m, amplitude of the opticalfield E = 10 V/m, and τ = 100 fs .First note that the magnitude of the current is of theorder of pA-nA which is within experimental reach. I x ismaximum when the polarization of light coincides withthe polar axis and decreases monotonically as the polar-ization turns away towards the y -axis. I y , on the otherhand, is nonmonotonic; it is zero when the light polariza-tion and the polar axis coincide, then rises to a maximumat 45 , and then decreases to zero for light polarized per-pendicular to the polar axis.In ultrafast pulsed experiments, the THz radiationemitted by the currents can be analyzed to study thenonlinear optical response of the system without need ofmechanical contact. In this scenario the system does nohave time to decay and the response is determined mainlyby the laser pulse characteristics not by the momentumdissipation mechanism. The above results indicate thatthe crystal structure, the geometry of the setup and lightpolarization can be used to uniquely characterize the jerkcurrent tensor components. Injection and shift currentshas been reported in THz spectroscopy in various mate-rials . XI. THIRD-ORDER INJECTION CURRENT
An explicit calculation of η is given in Appendix E.The result is η abcd (0 , ω, − ω,
0) = − πe (cid:126) V (cid:88) nm k f mn (cid:18) Ω adnm [ r bnm , r cmn ] + i [ r bnm , r cmn ; a ] ; d − iω nm ; a (cid:2) (cid:18) r dmn ω nm (cid:19) ; b , r cnm (cid:3)(cid:19) δ ( ω nm − ω ) − πe (cid:126) V (cid:88) nml k f mn ω nm ; a r dln ω nl [ r bnm , r cml ] D − ( ω nm , ω ) . (137)We defined Ω adnm ≡ Ω adn − Ω adm as the difference of Berryvector potentials. The Berry potential is related to theBerry curvature by Ω adn = (cid:80) e (cid:15) ade Ω en .The covariant derivative of r dmn /ω nm is with respectto the gauge dependent r dmn (see for example Eq. D6).The product r bnm r cmn ; a is gauge invariant and hence itscovariant derivative reduces to the standard derivative (see for example Eq. B1). To simplify notation we also defined[ O ( b ) , P ( c )] ≡ O ( b ) P ( c ) − O ( c ) P ( b ) , (138) D ± ( ω nm , ω ) ≡ δ ( ω nm − ω ) ± δ ( ω nm + ω ) , (139)where O, P are arbitrary matrix elements which dependon the cartesian indices b, c . For example[ r bnm , r cmn ] ≡ r bnm r cmn − r cnm r dmn . (140)One can see that η in Eq. 137 is manifestly anti-symmetric under exchange of b, c . In addition, it is6easy to show that η abcd (0 , ω, − ω,
0) is pure imaginaryand satisfies [ η abcd (0 , ω, − ω, ∗ = − η abcd (0 , ω, − ω,
0) = η abcd (0 , − ω, ω, b, c indicesimplies that η vanishes for linearly polarized light. η represents the current ddt J a (3)3 i = 6 (cid:88) bcd η abcd (0 , ω, − ω, E b ( ω ) E c ( − ω ) E d , (141)which varies as | J (3)3 i | ∼ η t, (142)in the absence of momentum relaxation and saturationeffects. A. Materials
In general, the 81 components of η are finite in bothcentrosymmetric and noncentrosymmetry crystal struc-tures. In practice, the symmetries of the 32 crystal classesgreatly reduce the number of independent components.For example, GaAs has ¯43 m point group, with 21 nonzerocomponents and 4 independent components. However, η is antisymmetric under exchange of b, c which reducesthe number of independent components to 1. In 2D ma-terials the number of components of η is also small. Forexample, single-layer GeS has mm η has only 2 independent components.In general, circular or unpolarized light will producethird-order injection current along the direction of thestatic field. Current transverse to the static field maynot be generated with unpolarized or linear polarization. B. Physical interpretation of third-order injectioncurrent
The presence of a static field gives rise to new physicalprocesses which we now describe in detail.
1. First term
The first term in Eq.(137) arises from the asymmetricinjection of carriers into anomalous velocity states. Tosee this, let us consider an electron wave packet in band n subject to a static field E d . The static field induces ananomalous contribution to the electron’s velocity whichgenerates a current given by (Eq. 56) J i, = − e (cid:126) V (cid:88) n k f n E × ΩΩΩ n , (143) - − i η ab c d ( A m / V s ) yyxxxxyy Photon energy (eV)
FIG. 5. Injection current response tensor η abcd in single-layerGeS near the band edges. η gives rise only to current trans-verse to the static field and vanishes for linearly polarizedlight. The tight-binding model parameters are described inSec. X F. where we used f n = ρ (0) nn . Taking a time derivative of theoccupations we obtain ddt J i, = − e (cid:126) V (cid:88) n k df n dt E × ΩΩΩ n . (144)This expression means that when the optical field isturned on electrons will be excited from the valence intoanomalous conduction states. To lowest order, i.e., sec-ond order in the optical field and first in the static field,Fermi’s Golden rule gives the one-photon injection rateshown in Eq.(89). Using Eqs.(89) we obtain ddt J a (3)3 i, = − πe (cid:126) (cid:88) b (cid:48) c (cid:48) d (cid:88) vc k Ω adcv r b (cid:48) cv r c (cid:48) vc δ ( ω cv − ω ) × E b (cid:48) ( ω ) E c (cid:48) ( − ω ) E d . (145)Using the fact that ω >
2. Second term
In the presence of a static field a wave packet driftsin the BZ giving rise to a current. Similarly, a pair ofwave packets could drift coherently in the presence of astatic field giving rise to a dipole current. To see this,consider the dipole velocity contribution to the currentin Eq. 94. Writing explicitly the small imaginary part ofthe external frequencies and taking the resonant part weobtain J a (2)3 i, = − iπe (cid:126) V (cid:88) bc (cid:88) nm k f nm (cid:2) r bnm ; a r cmn δ ( ω mn + ω )+ r cnm ; a r bmn δ ( ω mn − ω ) (cid:3) E b ( ω ) E c ( − ω ) . (146)7
03 2.0 2.5 3.00 I ( p A ) I ( p A ) Photon energy (eV) V Ixy V xy I (a)(b) yx x y FIG. 6. η -injection current in single-layer GeS near the bandedge. (a) shows the current parallel to the polar axis, I x and (b) the current transverse to the polar axis I y . Light iscircularly polarized and incident perpendicular to the planeof the GeS in both sample setups. Taking a time derivative of the dipole matrix elements,exchanging n, m indices, and making k → − k , we obtain ddt J a (3)3 i, = − iπe (cid:126) V (cid:88) bcd (cid:88) nm k f nm ∂∂k d (cid:0) r bnm ; a r cmn − r cnm ; a r bmn (cid:1) δ ( ω nm − ω ) E b ( ω ) E c ( − ω ) E d , (147)which can be recognized as the second term in Eq. 137.
3. Third term
The third term takes into account the change of theelectron distribution due to the static field. To see this,let us consider the current of an electron wave packet inband n to third order in the electric field. From Eq. 56 J a (3)3 i, = eV (cid:88) n k v an ρ (3) nn . (148)Taking a time derivative of the density matrix gives ddt J a (3)3 i, = eV (cid:88) n k v an ∂ρ (3) nn ∂t . (149)From Eq. 60 the time derivative of the density matrix is ∂ρ (3) nn ∂t = ei (cid:126) (cid:88) bm E b ( ρ (2) nm r bmn − r bnm ρ (2) mn ) − e (cid:126) (cid:88) b E b ρ (2) nn ; b . (150)Now consider the intraband part of the second order den-sity matrix obtained from Eq. 65 ρ (2) nm,i = ie (cid:126) (cid:88) bβcσ ¯ ρ (1) bβnm ; c ω mn − ω Σ E bβ E cσ e − iω Σ t , (151)where ω Σ = ω β + ω σ . The first order density matrix inthe presence of a static field is (see Eq. 63)¯ ρ (1) d nm = e (cid:126) f mn r dnm ω nm . (152)Substituting the above equations into Eq. 150 and takingthe resonant part we recover the third term in Eq. 137.The factor of two in Eq. 137 is due to two possible choicesfor the static electric field.
4. Fourth term
This contribution arises from electrons excited fromthe valence to conduction bands via an intermediate state l . These new states are generated by the presence ofstatic field. C. Third-order injection Hall current
Let us assume a static field is in the x -direction andcompute the current in the y -direction. An optical fieldof the form E = ˆ x E x ( ω ) e − iωt + ˆ y E y ( ω ) e − iωt + c.c. isincident perpendicular to the sample surface which wetake to define the xy -plane. From Eq. 141 the currenttransverse to the static field is ddt J y (3)3 iH = ς yx iH E x (153)where E a ( ω ) = | E a ( ω ) | e − iφ a and the Hall coefficient is ς yx iH ≡ iη yyxx | E x ( ω ) || E y ( ω ) | sin( φ x − φ y ) . (154)Similar to η , η vanishes for linear polarization φ x = φ y and is maximum for circularly polarized light. In a sim-ple relaxation time approximation, the dc singularity inthe conductivity (see Eq. 110) is cut off by a phenomeno-logical relaxation time τ as η − iω Σ → η τ − iω Σ . (155)8The current is J y (3)3 iH ∼ τ − iω Σ τ ς yx iH E x , (156)If ω Σ = 0, the current is proportional to τ the relax-ation of the diagonal elements of the density matrix. Forfrequencies larger than the Drude peak ω Σ τ (cid:29) D. Example: Third order injection current insingle-layer GeS
To get a sense of the η -injection current in real mate-rials we now calculate it for single-layer GeS. We use thesame 2-band, 2D tight-binding model of single-layer GeSand same sample geometry as in Sec. X F.Out of the 16 tensor components the antisymmetry inthe b, c indices and the mirror symmetry y → − y of thecrystal leaves only two independent components, yyxx and xxyy , shown in Fig. 5. These components allow cur-rent to flow only perpendicular to the static electric. Thesecond term in η is the dominant term followed by thethird, and the first terms which are one and two ordersof magnitude smaller respectively. We chose the optical field to be circularly polarizedand the static field is either along the polar axis of thesample or perpendicular to it E ( t ) = ˆ x E ( ω ) e − iωt + ˆ y E ( ω ) e − iωt + c.c., (157) E = ˆ x E x , or ˆ y E y , (158)where φ x − φ y = π/
2. The transverse currents are givenby I x (3)3 i = 12 Aτ iη xxyy | E x ( ω ) || E y ( ω ) | E x sin( φ y − φ x )(159) I y (3)3 i = 12 Aτ iη yyxx | E x ( ω ) || E y ( ω ) | E x sin( φ x − φ y )(160)where A = Ld is the transverse area of the sample. Notethat the current vanishes for linearly polarized light butis maximum for circular polarization. The calculated in-duced current is shown in Fig. 6a and 6b. Note that thephotocurrent is of the order of pA and of the same sign. XII. THIRD-ORDER SHIFT CURRENT
Explicit calculation of σ gives σ abcd (0 , ω, − ω,
0) = πe (cid:126) V (cid:88) nm k f mn (cid:2) { (cid:18) r dmn ω nm (cid:19) ; c , r bnm ; a } − { (cid:18) r dmn ; a ω nm (cid:19) ; c , r bnm } (cid:3) δ ( ω nm − ω ) − iπe (cid:126) V (cid:88) nml k f ln ω mn (cid:2) { r cnl , ( r dmn r blm ) ; a } − r dmn { r cnl ; a , r blm } (cid:3) D + ( ω nl , ω ) . (161)For details see Appendix F. In Eq. 161 we defined theanticummutator with respect to the b, c indices as { O ( b ) , P ( c ) } ≡ O ( b ) P ( c ) + O ( c ) P ( b ) (162)where O, P are arbitrary matrix elements. For example { (cid:18) r dmn ω nm (cid:19) ; c , r bnm ; a } ≡ (cid:18) r dmn ω nm (cid:19) ; c r bnm ; a + (cid:18) r dmn ω nm (cid:19) ; b r cnm ; a . (163) D + is defined in Eq. 139. Clearly, σ is symmet-ric under exchange of b, c , pure real, and satisfies σ abcd (0 , ω, − ω,
0) = σ acbd (0 , − ω, ω, J a (3) sh = 6 (cid:88) bcd σ abcd (0 , ω, − ω, E b ( ω ) E c ( − ω ) E d , (164) which, in the absence of momentum relaxation and satu-ration effects, is constant with illumination time (if quan-tum coherence is maintained). A. Materials
In general, the 81 components of σ are finite in bothcentrosymmetric and noncentrosymmetry crystal struc-tures. In practice, the symmetries of the 32 crystal classesgreatly reduce the number of independent components.For example, GaAs has ¯43 m point group, with 21 nonzerocomponents and 4 independent components . However, σ is symmetric under exchange of b, c which reduces thenumber of independent components to 3. In 2D mate-rials the number of components of σ is also small. Forexample, single-layer GeS has mm σ has only 6 independent components.In general, linear, circular or unpolarized light will pro-duce third-order shift current along the direction of thestatic field. Current transverse to the static field may notbe generated with unpolarized or circular polarization. B. Physical interpretation of the third-order shiftcurrent
1. First term
The first term in σ arises from the quantum inter-ference of the dipole velocity and interband band coher-ences. To see this, note that an oscillating external fieldcreates a dipole with wave packets in two distinct bands.Because the field oscillates the dipole velocity also oscil-lates, see Eq. 57. If the occupations of the bands, de-scribed by the density matrix, oscillate 180 out of phasewith respect to the velocity oscillations, a dc current canbe established. The process is mediated by the intrabandpart of the (second order) density matrix J a (3)3 sh, = − e (cid:126) V (cid:88) nm k E ( t ) · r nm ; a ρ (2) mn,i ( t ) , (165)where ρ (2) mn,i is the first term in Eq. 65 which clearly rep-resents the intraband part of ρ (2) mn . Setting one of thefields in ρ (2) to be static (say E dδ → E d ) we have J a (3)3 sh, = − ie (cid:126) V (cid:88) bβcσd (cid:88) nm k r bnm ; a ω mn − ω σ (cid:18) r dmn f nm ω mn (cid:19) ; c × E bβ E cσ E d e − iω Σ t (166)where ω Σ = ω β + ω σ . Symmetrizing with respect electricfield indices, substituting ω β = ± ω and ω σ = ∓ ω , andkeeping only resonant terms we recover the first term inEq. 161.
2. Second term
The second term in σ arises from the quantum inter-ference of interband coherences. To see this, note that astatic external field creates a dipole with wave packets intwo distinct bands. Because the field is static, the dipolevelocity is constant. The static dipole velocity togetherwith a nonoscillating interband occupation, can generatea dc current. This process is also mediated by the (static)intraband part of the (second order) density matrix J a (3)3 sh, = − e (cid:126) V (cid:88) d (cid:88) nm k E d r dnm ; a ρ (2) mn,i . (167) xxxxyyyyxyyxyxxyyyxx (x10) xxyy (x10) Photon energy (eV) σ ab c d ( x10 - A m / V ) FIG. 7. Shift current response tensor σ abcd of single-layer GeSnear the band edges. The model parameters are the same asin Sec. X F. The largest response is along the polarization axis x . The transverse response governed by xxyy and yyxx is anorder of magnitude smaller. Following the same procedure as above and after an in-tegration by parts it is easy to show that we recover thesecond term in Eq. 161.
3. Third and fourth term
The third and fourth terms in σ are not easily derivedfrom a simple model. These processes involve virtualtransitions to intermediate bands created by the staticexternal field and involve the interband part of the secondorder density matrix. C. Third-order shift Hall current
Let us assume a static field is in the x -direction andcompute the shift current in the y -direction. An opticalfield of the form E = ˆ x E x ( ω ) e − iωt + ˆ y E y ( ω ) e − iωt + c.c. is incident perpendicular to the sample surface which wetake as the xy -plane. The current transverse to the staticfield is J y (3)3 shH = ς yx shH E x (168)where E a ( ω ) = | E a ( ω ) | e − iφ a and the effective Hall con-ductivity is ς yx shH ≡ σ yyxx | E x ( ω ) || E y ( ω ) | cos( φ x − φ y ) . (169)Similar to σ , σ vanishes for circular polarization andis maximum for linear polarization φ x = φ y at 45 withrespect to the x -axis. Contrary to injection current, theshift current does not have a Drude-like dc divergencebut rather gives a finite contribution in this limit (whilequantum coherence is maintained).0 θ = 0 Photon energy (eV) I ( x10 - A ) I ( x10 - A ) V E E(t) θ xy θ = 0 yx I x I y
FIG. 8. σ -shift current in single-layer GeS near the band edgefor light linearly polarized at various angles θ with respect tothe polar axis. (a) current parallel to the polar axis I x islargest when the light’s polarization lies along the polar axis,and (b) current transverse to the polar axis I y is largest at θ = 45 . D. Example: third-order shift current insingle-layer GeS
To get a sense of the third-order shift current in realmaterials we now calculate it for single-layer GeS. We usethe same setup and tight-binding model of single-layerGeS as in Sec. X F.Because of the mirror symmetry y → − y of the model,only six tensor components are independent. As seen inFig. 7, the strongest is along the polar axis of magnitude ∼ × − Am/V . The component transverse to thestatic electric field σ yyxx (see Sec. XII C) is an order ofmagnitude smaller.The sample is rectangular of dimensions L × L andthickness d = 2 .
56 ˚A and is biased by an external batteryof voltage V as seen in Fig. 2c. For concreteness let usassume the optical field is incident perpendicularly to theplane of single-layer GeS as E ( t ) = ˆ x E x ( ω ) e − iωt + ˆ y E y ( ω ) e − iωt + c.c., (170) E = ˆ x E x . (171)The longitudinal and transverse currents are I x (3) sh = 6 A ( σ xxxx | E x ( ω ) | + σ xyyx | E y ( ω ) | ) E x (172) I y (3) sh = 6 Aσ yyxx | E x ( ω ) || E y ( ω ) | cos( φ x − φ y ) E x (173) where E x ( ω ) = E ( ω ) cos θe − iφ x , E y ( ω ) = E ( ω ) sin θe − iφ y , θ is the angle with the x -axis,and A = Ld is the transverse area of the sample. Notethat the current along the polar x -axis is independentof the polarization of light and hence, it will not vanisheven for unpolarized light. The transverse component ofthe current, on the other hand, vanishes for circularlypolarized (and unpolarized) light and is maximum forlinearly polarized light.We choose the optical field to be linearly polarized( φ x = φ y ) at an angle θ with the polar axis as shownin the inset to Fig. 8a. The figure shows the currentalong x and y -axis induced as a function of θ . We as-sumed the same semiconductor parameters as before,e.g., L = 100 µ m, V = 1V, E x = V /L = 10 V/m, am-plitude of the optical field E = 10 V/m, and τ = 100fs.First note that the magnitude of the currents is of theorder of pA-fA. I x is maximum when the polarization oflight coincides with the polar axis and decreases mono-tonically as the polarization turns away towards the y -axis. I y , on the other hand, is nonmonotonic: it is zerowhen the polarization and the polar axis coincide, thenrises to a maximum at 45 and then decreases to zeroagain for light polarized perpendicular to the polar axis. XIII. GENERALIZATIONSA. Snap current
By power counting it is easy to see that the leading di-vergence of χ is of order ω − , and that it occurs when allbut two of the external frequencies are zero. Proceedingas before we calculate the corresponding response tensor ς abcde (0 , ω, − ω, , ς abcde = 2 πe (cid:126) V (cid:88) nm k f mn (cid:2) ω nm ; ade r bnm r cmn + 3 ω nm ; ad ( r bnm r cmn ) ; e + ω nm ; a ( r bnm r cmn ) ; de (cid:3) δ ( ω nm − ω ) . (174)The tensor is symmetric in the b, c indices and representsa third derivative of the nonlinear current d J a (4) sp dt = 4! (cid:88) bcde ς abcde (0 , ω, − ω, , E b ( ω ) E c ( − ω ) E d E e , (175)where E d , E e represent static fields. By analogy with aparticle’s third derivative of its velocity we dub it snap current. The current grows as ∼ t with illuminationtime in the absence of momentum relaxation and satura-tion effects. Hence, it is proportional the third power ofthe relaxation time τ J a (4) sp ∼ τ (cid:88) bcde ς abcde E b ( ω ) E c ( − ω ) E d E e . (176)Note that we can think of the snap current as a secondorder photoconductivity. B. Higher-order singularities
One can show that the leading physical divergence of χ ni represents, in general, the n − − iω Σ ) χ i = ι + ( − iω Σ ) η + ( − iω Σ ) σ + ... (177)( − iω Σ ) χ i = ς + ( − iω Σ ) ι + ( − iω Σ ) η + ... (178)( − iω Σ ) χ i = κ + ( − iω Σ ) ς + ( − iω Σ ) ι + · · · (179)( − iω Σ ) χ i = (cid:36) + ( − iω Σ ) κ + ( − iω Σ ) ς + · · · . (180)...These higher-order analogs of the injection current arenamed by analogy with the time derivatives of a par-ticle’s velocity, e.g., jerk , snap , crackle , pop ,...,etc. anddenote them by, ι , ς , κ , (cid:36) ,.. respectively. Their physi-cal origin is similar to the injection current namely therate of carrier injection at current carrying states at time-reserved points in the BZ is asymmetric creating a polardistribution.An alternative formulation is the Laurent series for χ ni (or χ n since χ ne is regular or σ ( n ) ) as χ ni = ∞ (cid:88) l = − n a l z l (181)where z = − iω Σ and a l = 0 for frequencies less than thegap. The residues a − = η , a − = σ , a − = ι , etc., areformally given by a l = 12 πi (cid:73) | z | = ρ χ ni dzz l +1 , (182) ρ is the radius of convergence of the 1 /z series. In thesecalculations the limit ρ → (cid:15) → (but ω Σ = 0), the series starts from l > − n . XIV. EXPERIMENTAL SIGNATURES
In real materials, the measured current will be limitedby momentum relaxation mechanisms due to collisionswith other electrons, phonons, or impurities. For weak disordered insulators we expect the dc divergence of theconductivity in Eq. 110 will be cut off by a relaxationtime constant as σ (3) = ι ( τ − iω Σ ) + η τ − iω Σ + σ + · · · (183)assuming quantum coherence (time scale τ ) is main-tained. Calculation of τ requires a microscopic model ofmomentum relaxation which will be presented elsewhere.We have estimated the current of each contribution as-suming it can be detected separately. This is a challengein itself as is well documented in the literature. Here wepropose to use ultrafast THz spectroscopy together withthe symmetry of the crystal, the geometry of the setup,and the polarization of light to isolate these components.In ultrafast experiments, momentum relaxation plays aminor role (at least at short time scales) and the mag-nitude of the current is determined by the parametersof the lasers. For example the shift current magnitudefollows the envelope of the pulse . Recently, thesecond-order injection, shift or both currents have beenreported via THz radiation . In Table IIIwe present a summary of the jerk, injection and shiftHall-like responses of single-layer of GeS near the bandedge. As we can see, either the dependence on polariza-tion, the linearity of the static field, the order of magni-tude of the induced current, or the momentum relaxationtime scale can be used to distinguish them apart.
XV. CONCLUSIONS
The second-order injection and shift currents arearchetypical examples of nontrivial carrier dynamics ininsulators and semiconductors. In this paper we revisitedthe derivation of the second-order BPVE adding Fermisurface contributions to the theory and proposed a mi-croscopic interpretation of various BPVEs based on thecoherent motion of pairs of wave packets in the presenceof electric fields.We also studied the photoconductivity, i.e., a pho-tocurrent second-order in an optical and first order in astatic field, from the perspective of the third order elec-tric polarization susceptibility. Three new bulk photo-voltaic effects are found. We dub them jerk, third-orderinjection and third-order shift currents, respectively. Thejerk current and third-order injection currents can bethought of as a higher order versions of the standardsecond-order injection current and have essentially thesame microscopic origin, namely, the asymmetric rateof population of current-carrying states at time-reversedpoints in the BZ. The presence of the electric field, how-ever, gives rise the new contributions due to the anoma-lous and dipole velocity which are absent in the second-order injection current.The third-order shift current can be thought as ahigher order version of the second-order shift current.2
TABLE III. Summary of nonlinear Hall-like responses of single-layer GeS near the band edge. A static electric field is presentalong x which is taken to define the polar axis of GeS. In addition, an optical electric field is incident perpendicular to theplane of the sample which defines the xy -plane. The Hall current is in the y -axis. The sample geometry is shown in Fig. 2cand the details are in Sec. X F. I ≡ inversion symmetry, no I ≡ no inversion symmetry. I x is the current along the x -axis. ∗ Forcomparison, η - and σ -current is given for the same parameters. w.r.t. stands for ‘with respect to’.Current Momentum Dependence I vs. Hall current Hall current Hall current Sign of Hall current Ref. I y relaxation on E x no I dependence on vanishes for maximum for I x , I y magnitude ∗ Eq.polarization polarization polarization ι -jerk τ linear I, no I cos( φ x − φ y ) circular, linear at 45 +,- 10 − A 136linear E ( t ) (cid:107) x, y w.r.t. x -axis η -injection τ linear I, no I sin( φ x − φ y ) linear circular +,+ 10 − A 160 σ -shift τ linear I, no I cos( φ x − φ y ) circular, linear at 45 +,+ 10 − A 173linear E ( t ) (cid:107) x, y w.r.t. x -axis η -injection τ No no I 10 − A ∗ σ -shift τ No no I 10 − A ∗ It involves the coherent motion of pairs of wave pack-ets across the BZ. We showed that all photocurrents canbe understood using semiclassical wave packet dynamicsand showed that generalizations to higher order BPVEsare possible. Explicit expressions for the photocurrentsamenable for first-principles computations are given.
XVI. ACKNOWLEDGMENTS
We thank J.E. Sipe, R.A. Muniz, Y. Lin and C. Aversafor useful discussions. We acknowledge support fromDOE-NERSC Contract No. DE-AC02-05CH11231.
Appendix A: List of identities
Some definitions used in this paper are: v nn ( k ) = (cid:104) u n | v | u n (cid:105) ≡ v n ( k ) , (A1) f n ≡ f ( (cid:15) n ( k )) , (A2) f nm ≡ f n − f m , (A3) ξξξ nm ( k ) ≡ i (cid:104) u n |∇∇∇ k | u m (cid:105) , (A4) r nm ( k ) ≡ ξξξ nm ( k ) , ( m (cid:54) = n ) (A5) r nn ( k ) ≡ ω nm ≡ ω n − ω m . (A7)They describe velocity matrix elements (A1), Fermi dis-tribution (A2), Fermi function differences (A3), Berryconnection (A4) , off-diagonal (A5) and diagonal dipolematrix elements (A6), respectively, and frequency banddifferences (A7). u n is the periodic part of the Blochwave function (spinor index contracted). The covariantderivative of the dipole matrix elements is defined as r nm ; a ≡ (cid:20) ∂∂k a − i ( ξ ann − ξ amm ) (cid:21) r nm , (A8)or generally of any Bloch matrix element O nm as O nm ; a ≡ (cid:20) ∂∂k a − i ( ξ ann − ξ amm ) (cid:21) O nm . (A9)We also defined the commutator and anticommutatorwith respect to the Cartesian indices b, c as[ O ( b ) , K ( c )] ≡ O ( b ) K ( c ) − O ( c ) K ( b ) (A10) { O ( b ) , K ( c ) } ≡ O ( b ) K ( c ) + O ( c ) K ( b ) (A11)where O, K are any Bloch matrix elements. Some iden-tities used in this paper are: ω n ( − k ) = ω n ( k ) (A12) ω n ; a ( − k ) = − ω n ; a ( k ) = − ∂∂k a ω n ( k ) (A13) v nm ( − k ) = − v mn ( k ) = − [ v nm ( k )] ∗ (A14) r nm ( − k ) = r mn ( k ) = ( r mn ( − k )) ∗ (A15) r nm ; a ( − k ) = − r mn ; a ( k ) = − ( r nm ; a ( k )) ∗ (A16) ω nm ; a ( k ) = v an ( k ) − v am ( k ) = − ω nm ; a ( − k )= ω mn ; a ( − k ) (A17)ΩΩΩ n ( − k ) = − ΩΩΩ n ( k ) = − (ΩΩΩ n ( k )) ∗ . (A18)They arise from the hermicity of operators and the as-sumptions of time-reversal invariance of the ground state. (cid:126) ω n and ΩΩΩ n denote the band energy and Berry curvatureof band n .3 Appendix B: Derivation of η and σ from Taylorexpansion of χ To compute η abc (0 , ω, − ω ) and σ abc (0 , ω, − ω ) fromEq. 6, start from Eq. 72 and symmetrize ( − iω Σ ) χ i withrespect to pair-wise exchanges of electric fields indices b, β and c, σ . Then write explicitly the small imaginarypart of frequencies, ω β → ω β + i(cid:15) , ω σ → ω σ + i(cid:15) and let1 / ( x − i(cid:15) ) = 1 /x + iπδ ( x ). Next, set ω β = ω + n β ω Σ , ω σ = − ω + n σ ω Σ , 1 = n β + n σ , and Taylor expand real partsup to first order in ω Σ . It is easy to show that the non-resonant terms cancel and we obtain Eq. 77 and 80 asclaimed. In this calculation we used( r cnm r bmn ) ; a = r cnm ; a r bmn + r cnm r bmn ; a = ∂∂k a ( r cnm r bmn ) (B1) and some identities listed in Appendix A. Note that theexpression r cnm r bmn is gauge invariant and hence the co-variant derivative reduces to the standard crystal mo-mentum derivative. Appendix C: Expansion of χ i Using Eqs. 39, 56, 65, and 68 the third order suscep-tibility χ abcd ( − ω Σ , ω β , ω σ , ω ∆ ) can be written as χ = χ e + χ i where χ e C = − (cid:88) nm k r anm ω mn − ω Σ (cid:34) ω mn − ω ) (cid:18) r bmn f nm ω mn − ω β (cid:19) ; c (cid:35) ; d − i (cid:88) nmp k r anm ω mn − ω Σ (cid:34) ω mn − ω (cid:32) r bmp r cpn f pm ω mp − ω β − r cpm r bpn f np ω pn − ω β (cid:33)(cid:35) ; d − i (cid:88) nmp k r anm ω mn − ω Σ (cid:32) r bmp f pm ω mp − ω β (cid:33) ; c r dpn ω mp − ω − r bmp ω pn − ω (cid:32) r bpn f np ω pn − ω β (cid:33) ; c − (cid:88) nmpl k r anm ω mn − ω Σ (cid:34) r dmp ω pn − ω (cid:32) r bpl r clp f lp ω pl − ω β − r cpl r bln f nl ω lp − ω β (cid:33) − (cid:32) r bml r clp f lm ω ml − ω β − r cml r blp f pl ω lp − ω β (cid:33) r dpn ω mp − ω (cid:35) (C1) χ i C ≡ (cid:88) r =1 χ ir = 1 ω ω (cid:88) nm k ω nm ; a f mn (cid:18) r bnm r cmn ω nm − ω β (cid:19) ; d − ω (cid:88) nm k ω nm ; a r dmn ω nm − ω (cid:18) r bnm f mn ω nm − ω β (cid:19) ; c − iω (cid:88) nml k ω nm ; a r dmn ω nm − ω (cid:18) r bnl r clm f ln ω nl − ω β − r cnl r blm f ml ω lm − ω β (cid:19) − iω ω Σ (cid:88) nm k Ω adnm r bnm r cmn f mn ω nm − ω β + 1 ω Σ (cid:88) nm k r dmn ; a ω nm − ω (cid:18) r bnm f mn ω nm − ω β (cid:19) ; c + iω Σ (cid:88) nml k r dmn ; a ω nm − ω (cid:18) r bnl r clm f ln ω nl − ω β − r cnl r blm f ml ω lm − ω β (cid:19) . (C2)We defined C ≡ e / (cid:126) V , Ω adnm ≡ Ω adn − Ω adm , ω Σ ≡ ω β + ω σ + ω ∆ and ω ≡ ω β + ω σ . These expressions still needto be symmetrized with respect to pair-wise exchange ofelectric field indices ( b, β ), ( c, σ ), ( d, ∆). We note that, itis easier to calculate χ i from the intraband current J a (3) i rather than from P (3) i .Eq. C2 has a distinguishable structure, see Fig 9. Thefirst three terms in χ i are derived from the combination v an ρ (3) nn (Eq. 56). By analogy with χ i (Eq. 72), we wouldexpect these terms to be injection current-type of contri-4 χ i = χ i1 + χ i2 + χ i3 + χ i4 + χ i5 + χ i6 ι η σ Bloch vel Anomalous vel Dipole vel
FIG. 9. Origin of the first, second,..., contributions to the expressions for ι -jerk (Eq. 120), η -injection (Eq. 137) and σ -shift(Eq. 161) response tensors. Each of the 6 terms in the χ i originates from either the Bloch velocity (first three), anomalousvelocity (fourth) and the dipole velocity (5th and 6th). Due to the structure of the poles in χ i , the Bloch velocity and dipolevelocity contribute to multiple response functions. butions with one caveat; the first term has no analog in χ i since it is proportional to three powers of frequency, ω − ω − and is the most divergent at zero frequency. Thesecond and third terms, proportional to ω − , seem stan-dard injection coefficients similar to the first term in χ i .The fourth term is proportional to ( ω Σ ω ) − and arisesfrom the anomalous velocity ( E × ΩΩΩ) a ρ (2) nn . It is an in-jection current-type of coefficient. The fifth and sixthterms, proportional to ω − , originate from E · r nm ; a ρ (2) mn and hence are expected to be shift current-type of con-tributions.The goal in the next three sections (D, E and F) is tocalculate the coefficients ι , η , σ in the expansion( − iω Σ ) χ i = ι + ( − iω Σ ) η + ( − iω Σ ) σ + · · · (C3)To avoid cumbersome notation, we write the susceptibil-ities with the additional factors as( − iω Σ ) χ abcd i C → χ i . (C4)The strategy is to parametrize (the real part of) the ex-ternal frequencies as ω β = ω + n β ω Σ ,ω σ = − ω + n σ ω Σ ,ω ∆ = 0 , (C5)subject to n β + n σ = 1. Fig. 9 summarizes the result. Appendix D: Derivation of ι ι derives from χ i and χ i .
1. First term of ι Integrate by parts χ i and symmetrize it with respectto pair-wise exchange of electric field indices ( b, β ), ( c, σ ),( d, ∆) to obtain χ i, ≡ (cid:88) l =1 χ i, ,l = − iω Σ (cid:88) nm k ω nm ; ad f mn r bnm r cmn ( ω nm − ω β )( ω nm + ω σ ) − iω Σ (cid:88) nm k ω nm ; ac f mn r bnm r dmn ( ω nm − ω β )( ω nm + ω ∆ ) − iω Σ (cid:88) nm k ω nm ; ab f mn r dnm r cmn ( ω nm − ω ∆ )( ω nm + ω σ ) . (D1)The second and third terms will cancel against otherterms as we show later, but the first term will contributeto ι . By partial fractions and writing explicitly the imag-inary parts of the frequencies, the first term gives χ i, , = − iω Σ ω β + ω σ ) (cid:88) nm k ω nm ; ad f mn r bnm r cmn ( ω nm − ω β − i(cid:15) ) − iω Σ ω β + ω σ ) (cid:88) nm k ω nm ; ad f mn r cnm r bmn ( ω nm − ω σ − i(cid:15) ) . (D2)Using Eq. C5, 1 / ( x − i(cid:15) ) = 1 /x + iπδ ( x ), and expandingin powers of ω Σ to first order we obtain χ i, , = 2 π (cid:88) nm k ω nm ; ad f mn r bnm r cmn δ ( ω nm − ω ) − iω Σ (cid:88) nm k ω nm ; ad f mn r bnm r cmn ∂∂ω (cid:18) ω nm − ω (cid:19) (D3)The first term is independent of ω Σ and vanishes for fre-quencies smaller than the energy band gap. This is thefirst term of ι in Eq. 120. The second nonresonant termwill cancel against other terms.5
2. Second term of ι This contribution is obtained from χ i . To see this,let us symmetrize the second term in Eq. C2. After twointegration by parts we obtain χ i ≡ (cid:88) l =1 χ i ,l = iω Σ (cid:88) nm k ω nm ; ac r dmn r bnm f mn ( ω nm − ω β − ω σ )( ω nm − ω β ) + iω Σ (cid:88) nm k ω nm ; a r bnm f mn ω nm − ω β (cid:18) r dmn ω nm − ω β − ω σ (cid:19) ; c + iω Σ (cid:88) nm k ω nm ; ab r dmn r cnm f mn ( ω nm − ω β − ω σ )( ω nm − ω σ ) + iω Σ (cid:88) nm k ω nm ; a r cnm f mn ω nm − ω σ (cid:18) r dmn ω nm − ω β − ω σ (cid:19) ; b − iω Σ (cid:88) nm k ω nm ; a r cmn f mn ω nm − ω β − ω ∆ (cid:18) r bnm ω nm − ω β (cid:19) ; d − iω Σ (cid:88) nm k ω nm ; a r bmn f mn ω nm − ω ∆ − ω σ (cid:18) r dnm ω nm − ω ∆ (cid:19) ; c − iω Σ (cid:88) nm k ω nm ; a r bmn f mn ω nm − ω σ − ω ∆ (cid:18) r cnm ω nm − ω σ (cid:19) ; d − iω Σ (cid:88) nm k ω nm ; a r cmn f mn ω nm − ω ∆ − ω β (cid:18) r dnm ω nm − ω ∆ (cid:19) ; b . (D4)There are eight terms. To O ( ω Σ ), the l = 1 , l = 2 , l = 4 , η in Eq. 137 (see next section). The l = 5 , ι .Note that we can set ω β + ω σ = 0 (or ω ∆ = 0) (wherethis combination appears) since the pairs of poles in theseexpressions are distinct. This is not true in l = 5 , l = 5 term we obtain χ i , = − iω Σ (cid:88) nm k ω nm ; a r cmn r bnm ; d f mn ( ω nm − ω )( ω nm − ω β )+ iω Σ (cid:88) nm k ω nm ; a r cmn r bnm f mn ω nm ; d ( ω nm − ω )( ω nm − ω β ) (D5)here ω = ω β + ω ∆ and we used (cid:18) r dnm ω nm − ω ∆ (cid:19) ; c = r dnm ; c ω nm − ω ∆ − r dnm ω nm ; c ( ω nm − ω ∆ ) . (D6)Now obtain simple poles via partial fractions. The termwith a square of frequencies in denominator can be han-dled by ω nm ; d ( ω nm − ω ∆ ) = − ∂∂k d ( ω nm − ω β ) − , (D7)and a partial integration. Next, write the imaginary partof frequencies, use 1 / ( x − i(cid:15) ) = 1 /x + iπδ ( x ), and set ω β = ω + n β ω Σ , ω σ = − ω + n σ ω Σ , and 1 = n β + n σ . Note thatwith these definitions ω = ω + (1 + n β ) ω Σ . Now expandto second order in ω Σ and set (without expanding) ω ∆ = ω Σ . After some algebra we obtain χ i , = iω Σ (cid:88) nm k r cmn r bnm ; d f mn ∂∂k a (cid:18) ω nm − ω (cid:19) + iω Σ (cid:88) nm k ∂∂k d ( ω nm ; a r cmn r bnm ) f mn ∂∂ω (cid:18) ω nm − ω (cid:19) + π (cid:88) nm k ∂∂k d ( ω nm ; a r cmn r bnm ) f mn δ ( ω nm − ω ) (D8)In this calculation we have used the identity ∂∂ω (cid:18) ω nm − ω (cid:19) = − ∂∂ω nm (cid:18) ω nm − ω (cid:19) . (D9)Note that the third term in (D8) contributes to ι . Theother two nonresonant terms will eventually cancel. Asimilar calculation for the l = 7 terms gives χ i , = iω Σ (cid:88) nm k r bmn r cnm ; d f mn ∂∂k a (cid:18) ω nm + ω (cid:19) − iω Σ (cid:88) nm k ∂∂k d ( ω nm ; a r bmn r cnm ) f mn ∂∂ω (cid:18) ω nm + ω (cid:19) + π (cid:88) nm k ∂∂k d ( ω nm ; a r bmn r cnm ) f mn δ ( ω nm + ω ) (D10)Combining the l = 5 and l = 7 terms above and usingB1 we obtain χ i , + χ i , = iω Σ (cid:88) nm k ω nm ; ad r cmn r bnm f mn ∂∂ω (cid:18) ω nm − ω (cid:19) + 2 π (cid:88) nm k ∂∂k d ( ω nm ; a r cmn r bnm ) f mn δ ( ω nm − ω ) (D11)6The first term is nonresonant and will cancel against thesecond term in Eq. D3. The second term combined withthe first term in Eq. D3 gives ι in Eq. 120. Appendix E: Derivation of η We now derive each of the contributions to η inEq. 137.
1. First term of η The first term in η comes from χ i . Symmetrizing χ i in Eq. D4 and after partial fractions we obtain χ i, ≡ (cid:88) l =1 χ i, ,l (E1)= ω ω β + ω σ ) (cid:88) nm k Ω adnm f mn r bnm r cmn (cid:20) ω nm − ω β − ω nm + ω σ (cid:21) + ω (cid:88) nm k Ω acnm f mn r bnm r dmn ( ω nm − ω β )( ω nm + ω ∆ )+ ω (cid:88) nm k Ω abnm f mn r dnm r cmn ( ω nm − ω ∆ )( ω nm + ω σ ) (E2) Only the first term contributes to η . Writing the imag-inary parts of the frequencies, setting ω β = ω + n β ω Σ , ω σ = − ω + n σ ω Σ , and Taylor expanding, we obtain toleading order in ω Σ χ i, , = 2 iπω Σ (cid:88) nm k Ω adnm f mn r bnm r cmn δ ( ω nm − ω )+ ω (cid:88) nm k Ω adnm f mn r bnm r cmn ∂∂ω (cid:18) ω nm − ω (cid:19) (E3)Adding 1/2 of the first term to 1/2 of itself and letting k → − k in the second term we obtain the first contribu-tion of η in Eq. 137. The second term cancels againstother nonresonant contributions.
2. Second term of η This term arises from χ i . Symmetrizing we obtain χ i ≡ (cid:88) l =1 χ i ,l = iω (cid:88) nm k r dmn ; a f mn ω nm − ω β − ω σ (cid:18) r bnm ω nm − ω β (cid:19) ; c + iω (cid:88) nm k r dmn ; a f mn ω nm − ω β − ω σ (cid:18) r cnm ω nm − ω σ (cid:19) ; b + iω (cid:88) nm k r cmn ; a f mn ω nm − ω β − ω ∆ (cid:18) r bnm ω nm − ω β (cid:19) ; d + iω (cid:88) nm k r bmn ; a f mn ω nm − ω ∆ − ω σ (cid:18) r dnm ω nm − ω ∆ (cid:19) ; c + iω (cid:88) nm k r bmn ; a f mn ω nm − ω σ − ω ∆ (cid:18) r cnm ω nm − ω σ (cid:19) ; d + iω (cid:88) nm k r cmn ; a f mn ω nm − ω ∆ − ω δ (cid:18) r dnm ω nm − ω ∆ (cid:19) ; b . (E4)Let us consider χ i , first χ i , = iω (cid:88) nm k r cmn ; a f mn ω nm − ω (cid:18) r bnm ω nm − ω β (cid:19) ; d , (E5)where ω = ω β + ω ∆ . Performing a partial fraction ex-pansion, a substitution 1 / ( x − i(cid:15) ) = 1 /x + iπδ ( x ), fol-lowed by a Taylor expansion (to second order) in ω Σ of the real part about ( ω β , ω σ ) = ( ω, − ω ) using ω β = ω + n β ω Σ , ω σ = − ω + n σ ω Σ such that ω = ω +(1+ n β ) ω Σ ,we obtain χ i , = iω (cid:88) nm k r cmn ; a r bnm ; d f mn ∂∂ω (cid:18) ω nm − ω (cid:19) − iω (cid:88) nm k ( r cmn ; a r bnm ) ; d f mn ∂∂ω (cid:18) ω nm − ω (cid:19) + iω Σ ( iπ )6 (cid:88) nm k ( r cmn ; a r bnm ) ; d f mn δ ( ω nm − ω ) (E6)the first two terms are nonresonant contributions whichcancel against other terms. A similar analysis of χ i , gives7 χ i , = − iω (cid:88) nm k r bmn ; a r cnm ; d f mn ∂∂ω (cid:18) ω nm + ω (cid:19) + iω (cid:88) nm k ( r bmn ; a r cnm ) ; d f mn ∂∂ω (cid:18) ω nm + ω (cid:19) + iω Σ ( iπ )6 (cid:88) nm k ( r bmn ; a r cnm ) ; d f mn δ ( ω nm + ω ) (E7)the first two terms are nonresonant contributions whichcancel against other terms. After changing indices n, m and k → − k we see that the third term in Eq. E6 plusthe third term in Eq. E7 gives the second term of η inEq. 137.
3. Third term of η The third contribution to Eq. 137 arises from χ i , + χ i , + χ i , + χ i , in Eq. D4. Note that we can set ω β + ω σ = 0 from the outset since the poles in theseexpressions are distinct. Setting 1 / ( x − i(cid:15) ) = 1 /x + iπδ ( x ) and Taylor expanding about ( ω β , ω σ ) = ( ω, − ω ) we seethat to leading order the nonresonant parts vanish andwe obtain χ i , + χ i , = − ω Σ π (cid:88) bm k ω nm ; a (cid:18) r dmn ω nm (cid:19) ; c r bnm f mn δ ( ω nm − ω ) . (E8)Similar manipulations lead to vanishing nonresonantterms and to χ i , + χ i , = − ω Σ π (cid:88) bm k ω nm ; a (cid:18) r dmn ω nm (cid:19) ; b r cnm f mn δ ( ω nm + ω ) . (E9)Relabeling of indices n, m , setting k → − k , and addingto Eq. E8 we recover the third term of η .
4. Fourth term of η The fourth term arises from χ i . Let us label the 12terms obtained after symmetrization of χ i as χ i ≡ (cid:88) l χ i ,l = ω Σ (cid:88) nmo k ω nm ; a r dmn ω nm − ω β − ω σ (cid:20) r bno r com f on ω no − ω β − r cno r bom f mo ω om − ω β (cid:21) + ω Σ (cid:88) nmo k ω nm ; a r dmn ω nm − ω σ − ω β (cid:20) r cno r bom f on ω no − ω σ − r bno r com f mo ω om − ω σ (cid:21) + ω Σ (cid:88) nmo k ω nm ; a r cmn ω nm − ω β − ω ∆ (cid:20) r bno r dom f on ω no − ω β − r dno r bom f mo ω om − ω β (cid:21) + ω Σ (cid:88) nmo k ω nm ; a r bmn ω nm − ω ∆ − ω σ (cid:20) r dno r com f on ω no − ω ∆ − r cno r dom f mo ω om − ω ∆ (cid:21) + ω Σ (cid:88) nmo k ω nm ; a r bmn ω nm − ω σ − ω ∆ (cid:20) r cno r dom f on ω no − ω σ − r dno r com f mo ω om − ω σ (cid:21) + ω Σ (cid:88) nmo k ω nm ; a r cmn ω nm − ω ∆ − ω β (cid:20) r dno r bom f on ω no − ω ∆ − r bno r dom f mo ω om − ω ∆ (cid:21) . (E10)We analyze the structure of χ i by dividing its termsinto two groups. The first group composed of the l =1 , , , l = 5-12 terms. The l = 5 , , ,
10 terms have pairs ofpoles separable by partial fractions and can be combinedwith the l = 12 , , , ω Σ , it is useful to notewe can set ω β + ω σ = 0 or ω ∆ = 0 in all terms from theoutset. This is because the poles in each term are alwaysdistinct and separable by simple partial fractions. Thisshould be contrasted with the l = 5 , l = 1 , χ i , + χ i , == ω Σ (cid:88) nmo k ω nm ; a r dmn r bno r com f on ω nm F + ( ω no , ω β ) , (E11)where F + is defined as F + ( ω no , ω β ) ≡ ω no − ω β − i(cid:15) + 1 ω no + ω β + i(cid:15) = H + ( ω no , ω β ) + iπD − ( ω no , ω β ) , (E12)and8 H ± ( ω no , ω β ) ≡ ω no − ω β ± ω no + ω β D ± ( ω no , ω β ) ≡ δ ( ω no − ω β ) ± δ ( ω no + ω β ) . (E13)Similar manipulations for the sum of the l = 3 , χ i , + χ i , == ω Σ (cid:88) nmo k ω nm ; a r dmn r cno r bom f on ω nm F + ( ω no , ω σ ) . (E14)Adding the l = 1-4 contributions we find (cid:88) l χ i ,l == ω Σ (cid:88) nmo k ω nm ; a r dmn ω nm ( r bno r com + r cno r bom ) f on H + ( ω no , ω )+ iπω Σ (cid:88) nmo k ω nm ; a r dmn ω nm ( r bno r com − r cno r bom ) f on D − ( ω no , ω )(E15) The first term will cancel against other nonresonant con-tributions.Next we consider the group of l = 5 , , ,
10. It is easyto show these terms can be written as χ i , ≡ (cid:88) l χ i , ,l = ω Σ (cid:88) nmo k ω nm ; a r cmn r bno r dom f on ω mo (cid:20) ω nm − ω + iπδ ( ω nm − ω ) − ω no − ω − iπδ ( ω no − ω ) (cid:21) , (E16) χ i , ≡ (cid:88) l χ i , ,l = ω Σ (cid:88) nmo k ω nm ; a r cmn r dno r bom f mo ω no (cid:20) ω nm − ω + iπδ ( ω nm − ω ) − ω om − ω − iπδ ( ω om − ω ) (cid:21) , (E17) χ i , ≡ (cid:88) l χ i , ,l = ω Σ (cid:88) nmo k ω nm ; a r bmn r cno r dom f on ω mo (cid:20) ω nm + ω + iπδ ( ω nm + ω ) − ω no + ω − iπδ ( ω no + ω ) (cid:21) , (E18) χ i , ≡ (cid:88) l χ i , ,l = ω Σ (cid:88) nmo k ω nm ; a r bmn r dno r com f mo ω no (cid:20) ω nm + ω + iπδ ( ω nm + ω ) − ω om + ω − iπδ ( ω om + ω ) (cid:21) . (E19)We now combine them with the resonant ( r ) and non- resonant ( nr ) parts of the l = 12 , , , χ i , , + ( χ i , ) nr = − ω Σ (cid:88) nmo k ω nm ; a r cmn r bno r dom f mn ω om ( ω nm − ω ) , (E20) χ i , , + ( χ i , ) r = − iπω Σ (cid:88) nmo k ω nm ; a r cmn r bno r dom f mn ω om δ ( ω nm − ω ) , (E21) χ i , , + ( χ i , ) nr = ω Σ (cid:88) nmo k ω nm ; a r cmn r dno r bom f mn ω no ( ω nm − ω ) , (E22) χ i , , + ( χ i , ) r = − iπω Σ (cid:88) nmo k ω nm ; a r cmn r dno r bom f mn ω on δ ( ω nm − ω ) , (E23) χ i , , + ( χ i , ) nr = − ω Σ (cid:88) nmo k ω nm ; a r bmn r cno r dom f mn ω om ( ω nm + ω ) , (E24) χ i , , + ( χ i , ) r = − iπω Σ (cid:88) nmo k ω nm ; a r bmn r cno r dom f mn ω om δ ( ω nm + ω ) , (E25) χ i , , + ( χ i , ) nr = ω Σ (cid:88) nmo k ω nm ; a r bmn r dno r com f mn ω no ( ω nm + ω ) , (E26) χ i , , + ( χ i , ) r = iπω Σ (cid:88) nmo k ω nm ; a r bmn r dno r com f mn ω no δ ( ω nm + ω ) . (E27)Now we want to show that to O ( ω Σ ) the resonant partof the sum of the l = 1-4 and l = 5-12 groups givesthe fourth term of η and the nonresonant part vanishes.First the resonant contributions. a. Resonant contributions Let n ↔ m and k → − k in χ i , , and add to χ i , , to obtain χ i , , + χ i , , = − iπω Σ (cid:88) nmo k ω nm ; a r cmn r dom r bno f on ω mo D − ( ω no , ω ) . (E28)Similar manipulations on χ i , , and χ i , , give χ i , , + χ i , , = iπω Σ (cid:88) nmo k ω nm ; a r bmn r dom r cno f on ω mo D − ( ω no , ω ) . (E29)Adding Eq. E28 and E29 gives χ i , , + χ i , , + χ i , , + χ i , , = iπω Σ (cid:88) nmo k ω nm ; a r dom ω mo ( r bmn r cno − r cmn r bno ) f on D − ( ω no , ω ) . (E30) Performing analogous manipulations, add Eq. E21 toEq. E23 and Eq. E25 to Eq. E27 to obtain χ i , , + ( χ i , ) r + χ i , , + ( χ i , ) r = − iπω Σ (cid:88) nmo k ω nm ; a r cmn r bno r dom f mn ω om D − ( ω nm , ω ) , (E31)and χ i , , + ( χ i , ) r + χ i , , + ( χ i , ) r = iπω Σ (cid:88) nmo k ω nm ; a r bmn r dom r cno f mn ω om D − ( ω nm , ω ) . (E32)respectively. After n ↔ l , and k → − k in Eq. E30 addto Eq. E15 to obtain( (cid:88) l χ i ,l ) r + χ i , , + χ i , , + χ i , , + χ i , , = iπω Σ (cid:88) nmo k ω nl ; a r dmn ω nm ( r bno r com − r cno r bom ) f on D − ( ω no , ω ) . (E33)Now add Eq. E31 and Eq. E32 to obtain0 χ i , , + ( χ i , ) r + χ i , , + ( χ i , ) r + χ i , , + ( χ i , ) r + χ i , , + ( χ i , ) r = iπω Σ (cid:88) nmo k ω no ; a r dmn ω nm ( r bno r com − r cno r bom ) f on D − ( ω no , ω ) . (E34)Finally, the sum of all resonant terms in χ i to linearorder in ω Σ amounts to adding Eq. E33 to Eq. E34. Theresult is E
33 + E
34 =2 iπω Σ (cid:88) nmo k ω no ; a r dmn ω nm ( r bno r com − r cno r bon ) f on D − ( ω no , ω )(E35)which is the fourth term in η . b. Nonresonant contributions The sum of the (nonresonant) third terms in Eqs. E16,E17, E18 and E19 gives χ i , , + χ i , , + χ i , , + χ i , , = − ω Σ (cid:88) nmo k ω nm ; a r dom ( r cmn r bno + r bmn r cno ) f on ω mo H + ( ω no , ω ) . (E36)Next, the sum of Eqs. E20 and E22 and of Eq. E24 andE26 gives χ i , , + ( χ i , ) nr + χ i , , + ( χ i , ) nr = − ω Σ (cid:88) nmo k ω nm ; a r cmn r bno r dom f mn ω om H + ( ω nm , ω ) , (E37) χ i , , + ( χ i , ) nr + χ i , , + ( χ i , ) nr = − ω Σ (cid:88) nmo k ω nm ; a r bmn r cno r dom f mn ω om H + ( ω nm , ω ) . (E38)After l ↔ n and k → − k in Eq. E36 combined with thenonresonant part of Eq. E15 we obtain( (cid:88) l χ i ,l ) nr + χ i , , + χ i , , + χ i , , + χ i , , = ω Σ (cid:88) nmo k ω no ; a r dmn ( r com r bno + r bom r cno ) f on ω nm H + ( ω no , ω ) . (E39) Adding Eq. E37 and Eq. E38 we obtain E
37 + E
38 = − ω Σ (cid:88) nmo k ω nm ; a r dom ( r cno r bmn + r bno r cmn ) f mn ω om H + ( ω no , ω ) , (E40)which after l ↔ n and n ↔ m , is seen to cancel Eq. E39.This concludes the proof that to linear order on ω Σ thenonresonant terms vanish. Appendix F: Derivation of σ
1. First and second terms in σ Consider χ i , and χ i , in Eq.E4. In these terms wecan set ω β + ω σ = 0. Using 1 / ( x − i(cid:15) ) = 1 /x + iπδ ( x )and ∂∂k c (cid:32) r dmn ; a r bnm ω nm (cid:33) = (cid:32) r dmn ; a ω nm (cid:33) ; c r bnm + (cid:32) r dmn ; a ω nm (cid:33) r bnm ; c (F1)the resonant parts are( χ i , + χ i , ) r = πω (cid:88) nm k f mn (cid:2)(cid:0) r dmn ; a ω nm (cid:1) ; c r bnm + (cid:0) r dmn ; a ω nm (cid:1) ; b r cnm (cid:3) δ ( ω nm − ω )(F2)Similar manipulations on χ i , and χ i , in Eq. E4 yieldthe rest of the terms in the square brackets in σ . Thenonresonant parts can be shown to vanish.
2. Third and fourth terms in σ These contributions to σ arise from χ i in Eq. C2.It can be shown that the nonresonant parts vanish andthe resonant part gives the third and fourth term in σ .Since the algebraic steps are very similar to those usedin finding the third term in η we omit the derivation. Appendix G: two-band model of single-layer GeS
We consider a two-band, 2D model of single-layer GeSgiven by the Hamiltonian H = f σ + f a σ a , (G1)where σ a , a = x, y, z are the standard Pauli matrices and σ is the 2 × f a u c = A (cid:18) f x − if y (cid:15) − f z (cid:19) (G2) u v = A (cid:18) f z − (cid:15)f x + if y (cid:19) , (G3)where A − = 2 (cid:15) ( (cid:15) − f z ) is the normalization and eigen-values by E c,v = f ± (cid:15) where (cid:15) = √ f a f a and c, v denotethe conduction and valence band respectively. An arbi-trary phase factor in the eigenvectors has been omitted,since the final expressions are independent of this phase.The Bloch wave functions are constructed as ψ n k = (cid:88) R e i k · R [ u (1) n φ ( r − R )+ e i k · r u (2) n φ ( r − r − R )] , (G4)where u ( i ) n denotes the eigenvector corresponding toeigenvalue n = v, c (valence, conduction) and i = 1 , r = ( a , B with respect to site A which istaken to be the origin. φ ( r ) are p z -orbitals and R runs over all lattice positions. Notice that the phase of thewave function at site B is different than that at site A .The hopping parameters of the Hamiltonian are f = 2 t (cid:48) [cos k · a + cos k · a ]+ 2 t (cid:48) cos k · ( a − a ) , (G5) f x − if y = e i k · r ( t + t Φ k + t Φ ∗ k ) , (G6) f z = ∆ , (G7)where Φ k ≡ e − i k · a + e − i k · a , ∆ is the onsite potentialand t , t , t , t (cid:48) , t (cid:48) are hopping matrix elements as indi-cated in Fig. 2(c). a = ( a x , − a y ) , a = ( a x , a y ) are theprimitive lattice vectors.For single-layer GeS the parameters are: ( a x , a y , d ) =(4 . / , . / , .
56) ˚A, where d is the thickness ofthe slab, a = 0 .
62 ˚A, and ( t , t , t , t (cid:48) , t (cid:48) , ∆) =( − . , . , . , . , − . , .
41) eV. It was shownthat these parameters reproduce the band structure andgeometry of the wavefunction in the vicinity of theGamma point . To compare with bulk values the re-sults are multiplied by 2 /d . The factor of 2 takes intoaccount the smaller unit cell of the tight-binding model. R. Karplus and J. M. Luttinger, Phys. Rev. , 1154(1954). N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, andN. P. Ong, Rev. Mod. Phys. , 1539 (2010). N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev.Mod. Phys. , 015001 (2018). A. Vishwanath, Physics , 84 (2015). D. H. Auston, A. M. Glass, and A. A. Ballman, Phys.Rev. Lett. , 897 (1972). A. M. Glass, D. von der Linde, and T. J. Negran, AppliedPhysics Letters , 233 (1974). W. T. H. Koch, R. Munser, W. Ruppel, and P. Wrfel,Ferroelectrics , 305 (1976). V. I. Belinicher and B. I. Sturman, Physics-Uspekhi ,199 (1980). B. I. Sturman and V. M. Fridkin,
The Photovoltaic andPhotorefractive Effects in Non-CentrosymmetricMaterials (Gordon and Breach Science Publishers, Philadelphia,1992). V. I. Belinicher and B. I. Sturman, Ferroelectrics , 29(1988). R. von Baltz and W. Kraut, Phys. Rev. B , 5590 (1981). N. Laman, A. I. Shkrebtii, J. E. Sipe, and H. M. van Driel,Applied Physics Letters , 2581 (1999). N. Laman, M. Bieler, and H. M. van Driel, Journal ofApplied Physics , 103507 (2005). J. E. Sipe and A. I. Shkrebtii, Phys. Rev. B , 5337(2000). H. M. van Driel and J. E. Sipe, “Coherence control of pho-tocurrents in semiconductors,” (Springer, New York, NY,2001) Chap. 5, pp. 261–306. D. Cˆot´e, N. Laman, and H. M. van Driel, Applied PhysicsLetters , 905 (2002). A. Ghalgaoui, K. Reimann, M. Woerner, T. Elsaesser,C. Flytzanis, and K. Biermann, Phys. Rev. Lett. ,266602 (2018). M. Bieler, N. Laman, H. M. van Driel, and A. L. Smirl,Applied Physics Letters , 061102 (2005). M. Bieler, K. Pierz, U. Siegner, and P. Dawson, Phys.Rev. B , 161304 (2007). J. Rioux, G. Burkard, and J. E. Sipe, Phys. Rev. B ,195406 (2011). J. Rioux and J. Sipe, Physica E: Low-dimensional Systemsand Nanostructures , 1 (2012). C. Somma, K. Reimann, C. Flytzanis, T. Elsaesser, andM. Woerner, Phys. Rev. Lett. , 146602 (2014). M. Nakamura, F. Kagawa, T. Tanigaki, H. S. Park, T. Mat-suda, D. Shindo, Y. Tokura, and M. Kawasaki, Phys. Rev.Lett. , 156801 (2016). M. Holtz, C. Hauf, A.-A. Hern´andez Salvador, R. Costard,M. Woerner, and T. Elsaesser, Phys. Rev. B , 104302(2016). A. M. Rappe, I. Grinberg, and J. E. Spanier, Proceedingsof the National Academy of Sciences , 7191 (2017). J. E. Spanier, V. M. Fridkin, A. M. Rappe, A. R. Akba-shev, A. Polemi, Y. Qi, Z. Gu, S. M. Young, C. J. Hawley,D. Imbrenda, G. Xiao, A. L. Bennett-Jackson, and C. L.Johnson, Nature Photonics , 611 (2016). L. Z. Tan, F. Zheng, S. M. Young, F. Wang, S. Liu, andA. M. Rappe, npj Comput. Mater. , 16026 (2016). T. Rangel, B. M. Fregoso, B. S. Mendoza, T. Morimoto,J. E. Moore, and J. B. Neaton, Phys. Rev. Lett. , J. Iba˜nez Azpiroz, S. S. Tsirkin, and I. Souza, Phys. Rev.B , 245143 (2018). S. R. Panday, S. Barraza-Lopez, T. Rangel, andB. M. Fregoso, “Injection current in ferroelectric group-ivmonochalcogenide monolayers,” ArXiv:1811.06474 [cond-mat.mes-hall]. H. Wang and X. Qian, “Quantum nonlinear ferroic opticalhall effect,” ArXiv:1811.03133 [cond-mat.mes-hall]. B. M. Fregoso, T. Morimoto, and J. E. Moore, Phys. Rev.B , 075421 (2017). K. Kushnir, M. Wang, P. D. Fitzgerald, K. J. Koski, andL. V. Titova, ACS Energy Letters , 1429 (2017). M. Nakamura, S. Horiuchi, F. Kagawa, N. Ogawa, T. Ku-rumaji, Y. Tokura, and M. Kawasaki, Nature Communi-cation , 281 (2017). N. Ogawa, M. Sotome, Y. Kaneko, M. Ogino, andY. Tokura, Phys. Rev. B , 241203 (2017). K. Kushnir, Y. Qin, Y. Shen, G. Li, B. M. Fregoso, S. Ton-gay, and L. V. Titova, ACS Applied Materials & Interfaces , 5492 (2019). A. M. Burger, R. Agarwal, A. Aprelev, E. Schruba,A. Gutierrez-Perez, V. M. Fridkin, and J. E. Spanier, Sci-ence Adv. , eaau5588 (2019). M. Sotome, M. Nakamura, J. Fujioka, M. Ogino,Y. Kaneko, T. Morimoto, Y. Zhang, M. Kawasaki, N. Na-gaosa, Y. Tokura, and N. Ogawa, Applied Physics Letters , 151101 (2019). F. de Juan, A. G. Grushin, T. Morimoto, and J. E. Moore,Nature Communications , 15995 (2017). D. Rees, K. Manna, B. Lu, T. Morimoto, H. Borrmann,C. Felser, J. Moore, D. H. Torchinsky, and J. Orenstein,“Quantized photocurrents in the chiral multifold fermionsystem rhsi,” ArXiv:1902.03230 [cond-mat.mes-hall]. R. W. Boyd,
Nonlinear Optics (Academic Press; 2nd edi-tion, San Diego, USA, 2003). C. Aversa and J. E. Sipe, Phys. Rev. B , 14636 (1995). B. M. Fregoso, R. A. Muniz, and J. E. Sipe, Phys. Rev.Lett. , 176604 (2018). C. Aversa and J. E. Sipe, IEEE Journal of Quantum Elec-tronics , 1570 (1996). R. D. King-Smith and D. Vanderbilt, Phys. Rev. B ,1651 (1993). R. Resta, Rev. Mod. Phys. , 899 (1994). In the standard notation of susceptibilities a permittiv- ity of free space, (cid:15) , is factored out of χ n . For clarity ofnotation we dont factor this term. D. E. Aspnes, Phys. Rev. B , 4648 (1972). D. Culcer, A. Sekine, and A. H. MacDonald, Phys. Rev.B , 035106 (2017). M. Bass, P. A. Franken, J. F. Ward, and G. Weinreich,Phys. Rev. Lett. , 446 (1962). F. Nastos and J. E. Sipe, Phys. Rev. B , 235204 (2010). P. U. Jepsen, R. H. Jacobsen, and S. R. Keiding, J. Opt.Soc. Am. B , 2424 (1996). G. Li, K. Kushnir, M. Wang, Y. Dong, S. Chertopalov,A. M. Rao, V. N. Mochalin, R. Podila, K. Koski, and L. V.Titova, in (2018). E. I. Blount,
Solid State Physics: Advances in Researchand Applications , Vol. vol 13 (Academic Press, 1962). F. D. M. Haldane, Phys. Rev. Lett. , 206602 (2004). I. Sodemann and L. Fu, Phys. Rev. Lett. , 216806(2015). J. E. Moore and J. Orenstein, Phys. Rev. Lett. , 026805(2010). D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. ,1959 (2010). O. Matsyshyn and I. Sodemann, “The non-linear hallacceleration and the quantum rectification sum rule,”ArXiv:1907.02532 [cond-mat.mes-hall]. L. C. Gomes and A. Carvalho, Phys. Rev. B , 085406(2015). A. M. Cook, B. M. Fregoso, F. de Juan, S. Coh, and J. E.Moore, Nature Communications , 14176 (2017). S. R. Panday and B. M. Fregoso, Journal of Physics: Con-densed Matter , 43LT01 (2017). H. Wang and X. Qian, Nano Letters , 5027 (2017). D. Sun, C. Divin, J. Rioux, J. E. Sipe, C. Berger, W. A.de Heer, P. N. First, and T. B. Norris, Nano Letters ,1293 (2010). D. A. Bas, K. Vargas-Velez, S. Babakiray, T. A. Johnson,P. Borisov, T. D. Stanescu, D. Lederman, and A. D. Bris-tow, Applied Physics Letters , 041109 (2015). D. A. Bas, R. A. Muniz, S. Babakiray, D. Lederman, J. E.Sipe, and A. D. Bristow, Opt. Express , 23583 (2016). R. Atanasov, A. Hach´e, J. L. P. Hughes, H. M. van Driel,and J. E. Sipe, Phys. Rev. Lett.76