Consequences of Time-reversal-symmetry Breaking in the Light-Matter Interaction: Berry Curvature, Quantum Metric and Diabatic Motion

Nonlinear optical response is well studied in the context of semiconductors and has gained a renaissance in studies of topological materials in the recent decade. So far it mainly deals with non-magnetic materials and it is believed to root in the Berry curvature of the material band structure. In this work, we revisit the general formalism for the second-order optical response and focus on the consequences of the time-reversal-symmetry ($\mathcal{T}$) breaking, by a diagrammatic approach. We have identified three physical mechanisms to generate a dc photocurrent, i.e. the Berry curvature, the quantum metric, and the diabatic motion. All three effects appear in general for broken time-reversal symmetry and can be understood intuitively from the anomalous acceleration. The first two terms are respectively the antisymmetric and symmetric parts of the quantum geometric tensor. The last term is due to the dynamical antilocalization that appears from the phase accumulation between time-reversed fermion loops. Additionally, we derive the semiclassical conductivity that includes both intra- and inter-band effects. We find that $\mathcal{T}$-breaking generally enhances the conductivity by contributing the leading-order term ($\omega^{-2}$, where $\omega$ is the light frequency) while preserving $\mathcal{T}$ restricts the response to the term of next-to-leading order ($\omega^{-1}$).


I. INTRODUCTION
The Bulk photovoltaic effect (BPVE) [1][2][3] refers to the generation of a dc current from a uniform material by irridating strong light.Since the early 1980s, its understanding has been gradually established as the secondorder optical response [2,[4][5][6][7][8] that is closely related to the Berry phase [9], an intrinsic quantity of the material band structure.Recently it gained renewed interest in topological Weyl semimetals (WSMs) [10,11], where the Berry phase is believed to generate the giant photocurrent and second harmonics as well [12][13][14][15][16][17][18][19][20].Thus far, research on BPVE focuses on nonmagnetic materials.However, recent theoretical [21] and experimental [22][23][24] works on magnetic systems reveal distinct photocurrent and second harmonic generation, which cannot be merely derived from the Berry phase.Therefore, we are motivated to reexamine the second-order response theory and investigate the effects of time-reversal-symmetry (T ) breaking.
The non-linear current j created at the second order in the incident electric field E of light is defined as A focal point will be in the following the dc-current j(0; ω, −ω) which is created in response to irradiation of light.In a clean, gapped system with T , the BPVE is usually phrased in the form of two phenomena, respectively termed shift current and injection current.Both effects require inversion symmetry breaking, otherwise the momentum-space integral of the response function vanishes.The shift current denotes the current that is proportional to the positional shift (shift vector expressed in the form of the Berry connection) of the electron charge in an interband process.The shift current has previously been employed to successfully describe the BPVE in many compounds [2,6,[25][26][27][28].Injection current is a stronger response, which is described by a constant current creation with rate dj/ dt [5,[29][30][31][32][33][34][35].In systems with finite relaxation rates, dissipative processes will limit the current injection, thus rendering the nonlinear conductivity finite, with a magnitude proportional to the lifetime of the quasiparticles.In materials with T , the injection current will occur only in response to circular polarized light and a linear polarization can only lead to a shift current.This distinction is lost in compounds which break T , as pointed out by Ref. 21 very recently, where both shift and injection current might arise for either linear or circular polarized light.
Compared to the shift current generation, the intuitive understanding of current injection is less clear.One of the initial objectives of this work was to understand its physical origin.An important step was the realization that a single cone in a two-band Weyl model can produce injection current at a quantized rate, which is proportional to the topological charge of a Weyl cone [19].However, it seems unlikely for this identification of the injection current with topological charge to translate to the generic case.Also, an analogous procedure to express the shift current with the Berry curvature has not been reported.Given these findings, it is presently unclear to which extent the presence of Berry monopoles located at the Wey nodes determine the BPVE in WSMs.Indeed, topologically trivial semimetals like Graphene have zero Berry curvature and are still known to possess a whole assortment of unique transport properties, among them an abnormally large third order conductivity, merely due to their nodal points in the bandstructure and the surrounding linear Dirac cone dispersion [36][37][38].A positive statement about the topological origin of the large shift current observed for example in TaAS [12] thus requires an identification of this response with an observable related to the topological nature of the material.
Recently, an intriguing nonlinear anomalous Hall ef-arXiv:1911.05667v1[cond-mat.mes-hall]13 Nov 2019 fect [39][40][41][42] was derived in the semiclassical limit, which refers to intraband transitions in the small frequency case of the BPVE.It is attributed to the existence of the Berry curvature dipole in the semiclassical description of the second-order response [16,41] and shortly after discovered in the WSM, WTe 2 , layers [43][44][45][46].It is, however, unclear how the mechanisms for shift and injection currents, which are from interband transitions, are related to the semiclassical motion in terms of the Berry curvature dipole.
The present discussion also touches upon a practical question of how to properly model nonlinear response in realistic quantum systems with impurities or at finite temperature.Historically, a large part of the discussions are framed in the context of two gauge choices, the length gauge [5,6], and the velocity gauge [8,47].These names are derived from the way the applied electric field is coupled to the Hamiltonian.In the length gauge, electric field enters via a dipole energy which induces a polarization, while in the velocity gauge it is included through the electromagnetic gauge potential via minimal substitution in the momentum.Both approaches were originally conceived as complementary descriptions, but in the clean limit they yield identical results and can be related by a time-dependent gauge transformation [7,48].The situation is less clear-cut in realistic calculations using finite lifetimes [47,49].The issue is that non-linear optical response fundamentally involves interband processes which are difficult to model semiclassically.As such, there is some ambiguity with regards to the lifetimes which enter the expressions.While we will not attempt to solve this issue comprehensively, we offer a way to discuss both the adiabatic and the diabatic motion which arises under optical driving in terms of the so-called quantum geometric tensor [50].Since this object has a clear physical interpretation, it is possible to retrace which finite relaxation rates enter the expressions.
In the following, we proceed to systematically explore the second-order optical response where not only inversion symmetry is absent but also the T is broken.In particular, we derive formulas for the semiclassical response, shift current and injection current in terms of observable, physically distinct processes, some of which lead to previously unknown types of photocurrents as long as T is absent.Our starting point is the treatment of the BPVE in the diagrammatic approach [8] which can also account for finite lifetimes and offers a more immediate interpretation of the different contributions to second order response.Still, these expressions remain opaque to a semiclassical interpretation of the physical processes as described by the perturbation theory.Therefore, we propose to understand the nonlinear response by rewriting the second-order conductivity in terms of semiclassical accelerations, separately considering instantaneous and sequential multi-phonon processes, as dictated by the diagrammatics.
As the main results, we demonstrate that the secondorder optical response can be understood in terms of three distinct physical processes, each of them a part of the matrix elements of the anomalous quasiparticle acceleration in a Bloch band.While it is well known that the velocity matrix elements in band eigenbasis are not simply derivatives of the eigenenergies but contain an anomalous velocity originating from the nonzero Berry curvature [51], the anomalous acceleration has not yet been discussed at this same level of detail [26,[52][53][54].To be precise, consider a dispersion ε m with band index m, its momentum derivative v a mn for the spatial direction a and the Berry connection r a mn .The semiclassical center of mass motion of a wavepacket with position r and momentum k moving in band m subject to an electric field E has the equation of motion where the linear-response form of the anomalous velocity e E b Ω ba mm is due to the Berry curvature Ω ba mm = ∂ b r a mm − ∂ a r b mm .As a reminder, this result is most easily derived in the adiabatic approach by employing the identification ∂ t = (−e/ )E a ∂ ka , which follows from the semiclassical Lagrangian in minimal coupling [9].This means that the acceleration can be found by considering the expectation value ra However, this derivation of the acceleration is problematic if bilinear couplings to the electric field are to be retained, as this requires the equations for the semiclassical motion to be correct up to second order.Indeed, Gao et al. [55] reported some time ago that to second order in applied fields, there appears not only a polarization dependent renormalization of the dispersion ε(E), but also a field dependent positional shift a (E) of the semiclassical wave packet, both of which enter the equation of motion for a wavepacket, However, the discussion at the time was mostly concerned with the anomalous Hall effect and did not resolve these modifications with respect to second order optical response.We remedy this situation, but for the discussion at hand it turns out to be better to directly consider the matrix elements of the second momentum derivative Here, the last term is the normal derivative (i.e. a Drude term), while the other three give rise to the anomalous acceleration.These features of the anomalous motion appear in the second order optical response due to the close connection between the electric field and the derivative in time.However, a straightforward translation between anomalous motion and optical response is complicated by the fact that the perturbative response formulas intrinsically contain interband processes and thus off-diagonal matrix elements in the band basis.For this reason it went so far unnoticed that second-order response actually contains three types of anomalous contributions, unlike the linear response, where only the Berry curvature makes an appearance.The phenomenology is summarized in Tab.I. (i) The first term in Eq (4) is the well-studied adiabatic acceleration of the electron wavefunction due to the dipole interaction with the electric field.In the second-order response, this can always be formulated as a properly symmetrized triple product of velocity matrix elements, and it corresponds to the Berry curvature dipole term in the simplified semiclassical treatment.Importantly, this adiabatic acceleration appears in the shift current.
(ii) The second term is a diabatic effect, which can also be written as a velocity difference between occupied and unoccupied bands.This response type can be identified as a mismatch between the derivative ∂ kj v i of the quasiparticle velocity v i and the derivative ε∂ kj r i of the Berry connection r i .Speaking in terms of Feynman diagrams, the phases accumulated on momentum reversed (i.e.time reversed) paths do not cancel with each other.This phenomenon is responsible for current injection, which can therefore be viewed as the result of dynamical, i.e. finitefrequency antilocalization between time-reversed loops with three legs (see Fig. 1).
(iii) The third term contains what is known as the bare quantum metric.Somewhat surprisingly, it is impossible to write higher order conductivities only as a function of velocity matrix elements and the Berry curvature.This is manifest in the third type of anomalous acceleration, which is constituted by the symmetric part of the second derivative of the wave function, which is closely related but not identical to the gaugeinvariant quantum metric g ab [50,[56][57][58][59][60].In essence, second order response contains information about both the symmetric and antisymmetric part of the quantum geometric tensor (QGT), unlike the linear response coefficients which receive a contribution from the Berry curvature Ω ab mm only.We emphasize that both the Berry curvature and the bare quantum metric enter the second-order conductivity not directly, but in a form resembling a momentum derivative of these quantities.This is also the reason why it is necessary to consider the off-diagonal matrix elements of the acceleration, and why it is helpful to intermittently consider gauge-covariant and not exclusively gauge-invariant quantities.This nonwithstanding, our final expressions for the different contributions to the non-linear conductivity are gauge-invariant.(iv) Finally, we derive the semiclassical conductivity in the absence of T by employing a canonical expansion of the conductivity in powers of the applied frequency (ω).The leading order term (ω −2 ) appears only when breaking T and contains the above four contributions, as summarized in Table I.We find that nonlinear response contains terms of geometric origin exclusively in the next-leading order in frequency (ω −1 ).As T is broken, this next-leading order term is not given only by the Berry curvature dipole, but contains an additional diabatic term due to the intrinsically dissipative nature of optical transitions.
In addition, we report a robust way to obtain explicit sum rules for a general tight-binding Hamiltonian required to translate between velocity gauge and length gauge, even if the system is imperfect and quasiparticle lifetimes are finite.Our derivation shows that the sum rules should be viewed as an expansion of higher order vertices in terms of velocities and band energies.This is in agreement with the general formulation of sum rules in clean systems, where they appear in terms of vanishing nested covariant derivatives of an observable, as mandated by gauge invariance [7].However, the latter formulation is entirely formal, leaving the resolution of the derivatives of matrix elements untouched.It is precisely this point where additional lifetime dependencies appear, thus making the sum rules one possible source of disagreement between the length and the velocity gauge formalism.In particular, the often quoted resolution of the non-Abelian Berry connection in terms of band structure parameters as r a mn = v a mn /i mn does not hold if the quasiparticle lifetime is finite.In the past, another point of confusion was the use of commutator identities which are only valid for an infinite number of bands [7,61].For example, the relation generally does not hold for a tight-binding Hamiltonian.However, this problem is unrelated to gauging, it originates from the reliance on the relation p a = v a m, which is only true for an isotropic dispersion.The diagrammatic formalism does not necessitate assumptions of this kind, therefore the sum rules presented here do not suffer from these limitations.We find that the sum rules are only effective at reducing the zero frequency singularity of the second-order conductivity if the material is time-reversal symmetric.
The remaining sections are organized as follows.In Sec.II, we introduce the general formalism for nonlinear optical response for a clean system, relying mostly on the formulation in the velocity gauge.The introduction of finite lifetimes in a disordered system can be intuitively implemented within the diagrammatic approach.We discuss the phenomenology of this approach and also compare to the equivalent procedure in the length gauge.In Sec.III we explain the various processes which contribute to the second order response.Some time is spend on the semiclassical picture associated with current injection, and we derive an intuitive formulation of this phenomenon in terms of a phase difference accumulated on a closed path.We calculate the semiclassical conductivity including all multiband contributions in the absence of T .In Sec.IV, we illustrate our findings by numerical calculations on a concrete example of a WSM and a Dirac semimetal.This serves as a model for a system which has zero (Abelian) Berry curvature and yet exhibits a large BPVE.

II. PERTURBATION THEORY
We briefly recollect the basic ingredients needed to capture optical response when going beyond linear approximation.In the literature, this was achieved either by calculating the electric susceptibility or more directly by evaluating the current-current correlator to second order in the electric fields.These approaches are also known under the terms length gauge and velocity gauge.

A. Adiabatic formalism
The current created by the BPVE originates from the intraband motion of Bloch states in response to the irradiation.The landmark result for this effect was derived in [5] and reads Here, m, n are band indices, ε mn = ε m − ε n is the difference of the eigenvalues of the unpertubed Hamiltonian H 0 in units of frequency, i.e. ε m (k) = −1 m|H 0 |m ; the difference of Fermi factors is f mn = f m − f n and finally there is the difference of velocities, ∆ a mn = v a mm − v a nn , with v a mn = m|∂ ka H 0 |n .The non-Abelian Berry connection is r a mn = i m|∂ ka n with the generalized derivative r b nm,c = ∂ kc r b nm −i(r c nn −r c mm )r b nm , where r c nn = i n|∂ kc n is the Abelian Berry connection.The integration is over the first Brillouin zone, with k ≡ F BZ d 3 k/(2π) 3 .The current created through the bulk photovoltaic effect formally corresponds to an infinite polarizability, i.e. charge imbalances are continuously created.However, the first term in Eq. ( 8) for the current is infinite by itself, meaning that the second-order response may even produce an infinite current, a phenomenon termed current injection, i.e. current is continuously created.Originally, Eq. ( 8) was derived from a painstaking second order calculation of the adiabatic time evolution of the instantaneous eigenstates, also known as a reduced density matrix formalism.While this approach is conceptually very clear-cut, it relies on two important factors.On the one hand, the optical processes must be fast enough for lifetime effects to be negligible.This assumption becomes questionable for optical transitions involving resonances with transient states in highly excited levels.In the present context, one would expect the current injection to strongly depend on lifetime effects.Secondly, one needs to appropriately resolve the derivative r b nm,c by a sum rule.The latter task was recently reexamined in [27], finding for n = m Here, the second derivative w ab mn = m|∂ ka ∂ k b H 0 |n makes an appearance.This term is was originally not included because a quadratic Hamiltonian implies w ab nm = 0 for n = m.Generically w is only one of several higher order vertices which appear in an order-by-order analysis [8].
Properly including lifetime effects necessitates a closer look at the theoretical underpinnings of the adiabatic formalism.The starting point is a formulation of the problem in the length gauge.Equivalently, it follows from the proposition that the response of the crystal can be approximated by a dipole Hamiltonian [5] where the renormalized dispersion εnm (k, t) = ε n δ nm − e −1 r a nm E a (t) contains the Berry connection r.While this construction encapsulates the polarization response in a clean system, it relies on the identification of the adiabatic intraband motion with the motion induced by the electromagnetic field [9].Namely, the crystal momentum q in the isolated system and the gauge invariant momentum k are by definition related by k = q + eA(t)/ .This implies k = −eE(t)/ under the condition that q = 0, i.e. if the crystal momentum is conserved.

B. Diagrammatic formulation
In the velocity gauge, one instead considers the coupling to the gauge field through minimal substitution of the momentum operator p → p − eA/ in the unperturbed Hamiltonian H 0 .Expanding in powers of A yields therefore Here, we wrote everything in terms of operators and not their matrix elements, which would result in a very different expression due to the nature of the (Berry) covariant derivative in k-space, which is defined through the connection where O is a Hermitian operator.This connection plays a major role in the calculation of Eq. ( 11).The current operator follows from a variation with respect to the electromagnetic field, Henceforth we will refer to derivatives in k-space by ∂ a = ∂/∂ ka .The various matrix elements are therefore v a mn = m|∂ a H 0 |n , w ab mn = m|∂ a ∂ b H 0 |n and m|∂ a ∂ b ∂ c H 0 |n = u abc mn .We note that using the connection, Eq. ( 12), the velocity operator v a = ∂ a H 0 can be expanded in the eigenbasis to yield the familiar expressions for the anomalous velocity including the Berry curvature term.The same holds for all higher order derivatives.We will defer this step until the discussion of the semiclassical motion.
In total, in second order the conductivity receives contributions from four types of diagrams, two of which have two variants.They are depicted in Fig. 1.Diagrams (I) and (II) corresponds to the type of response originally considered in the earliest works [2,4].Importantly, they only carry velocity matrix elements.Diagrams (III)-(VI) contain higher derivatives in k, which in the original formulation of Eq. ( 8) appear through ∂ c r b nm , via sum rules.In the length gauge these higher order derivatives follow from an expansion of the time evolution of the instantaneous eigenstates beyond linear approximation [5,9].
The complete non-linear conductivity for clean system reads where ω = ω 1 + ω 2 .We point out that Eq. ( 14) differs from Ref. [8] by a sign for the band energy differences in diagrams (III)-(V) [62].As it is visible from the diagrams, the vertices w and u describe instantaneous processes, where two or three interaction events coincide.In contradistinction, diagrams (I) and (II) describe sequential interactions.We repeat that Eq. ( 8) was computed from the dipole Hamiltonian and necessarily assumes the quasiparticle motion to be uninterrupted.Thus it cannot differentiate between instantaneous and successive twophoton processes.In the diagrammatics, these processes are easily distinguishable from their vertex structure.This distinction can be useful in particular to tell apart effects of the anomalous velocity from the ones produced by higher derivatives.

C. Finite relaxation rates
It is possible to prove the equivalence of length and velocity gauge in the adiabatic limit (γ → 0) [7].On a formal level, the equivalence only requires a proper implementation of the time-dependent unitary transformation which translates between state vectors in scalar potential gauge (length gauge) and velocity gauge.
On the other hand, it is not obvious how to include finite relaxation rates into nonlinear response formulas.Dissimilar schemes have been proposed for length and velocity gauge, and in some cases these small changes in the inserted relaxation rates have led to markedly different results for the non-linear response.One possibility is to retain different relaxation rates for the non-equilibrium response in each order of the electromagnetic field.Generically, this will produce divergences at resonances due to a mismatch in rate differences [49,63].A careful analysis of this issue was performed recently in Ref. [47].In the reduced density matrix formalism a relaxation rate γ is introduced through the appended equation of motion for the density matrix ρ, where ρ 0 is the equilibrium state, Translated into velocity gauge this entails a relaxation of the system towards a modified equilibrium distribution ρ(A) [7].It is not obvious to which physical situation such a relaxation mechanism corresponds to.On the other hand, starting from diagrammatics, a commonly employed device is adiabatic switching, which regularizes the poles in the energy denominators by the replacement ω → ω + iγ.Most importantly, this means inserting 2γ (3γ, . . . ) in the formulas for second, (third, . . . ) harmonic generation whenever 2ω (3ω, . . . ) appears.Backporting this prescription to the length gauge, the authors in [47] present a case where a highly spiked resonance peak for third harmonic generation becomes noticeably regularized by just this replacement.
Why is it so difficult to establish a unique procedure?To demonstrate the difficulties associated with the standard treatment of second order response, let us try to interpret the following second order term written in the length gauge Taking ω = 0, this can be shown to assume the form of the Berry dipole, ∂ c Ω ab mm [41].However, performing the momentum derivative first, one obtains instead Looking at this, it remains unclear which physical process is at play, and the reduction to the Berry dipole by any means other than reverting the derivative is contrived.These issues hold for all such terms, the presence of derivatives in the length gauge formulation make it hard to do power counting in ω and it is equally laborious to keep track of the proper gauge invariant combinations of terms [5].The same issues persist for any finite imaginary parts.
Fundamental to the problem is the resolution of the velocity matrix elements in the eigenbasis, the defining relation which allows to translate response formulas between length and velocity gauge.In the adiabatic limit, one employs r mn = v a mn /i mn (m = n).The origin of this is the relation which is invariably true -it follows from the mere fact that the covariant derivative defines a connection in momentum space [8,9].However, for a small but nonzero relaxation rate Γ mn for excited states, the commutator in this expression is only approximately diagonal in the eigenbasis, and yields Note that Γ mn is symmetric.The off-diagonal velocity matrix elements are no longer just proportional to the Berry connection.Only in a simplified picture where the relaxation is happening purely between bands m and n, i.e.where the relaxation is intraband only, it holds for m = n that r a mn = v a mn /i[ mn + i(Γ mm − Γ nn )].For the reasons detailed above, we employ the orthodox view of renormalized perturbation theory to approach this issue in the diagrammatic formulation with velocity gauge.Namely, we analyze when the propagators are on-shell or off-shell in the respective diagrams, and use this to introduce relaxation rates according to the self-energies associated with low-energy or high-energy bands.The approach is physically transparent and does not rely on artificially setting all relaxation rates for all states to be the same.It also sidesteps the issues associated with quasiparticle decay in a density matrix formulation, instead relying on the Green's function formalism to implement finite lifetimes.More importantly, different relaxation rates appear depending on the type of process and not on the order in perturbation theory.This also means that second harmonic generation does not involve the same relaxation rates as dc-current generation.For example, the three-point loop depicted in Fig 2 involves two consecutive photon processes, which means the quasiparticle partially propagates in highly excited bands.However, comparing the on-shell processes (Ω = 0), for the photovoltaic effect (ω 1 = −ω 2 ) the quasiparticle propagators are dissimilar to the case of second harmonic generation (ω 1 = ω 2 ).It seems consequent to take this into account by using appropriate (i.e.dissimilar) relaxation rates.In the calculation, this reasoning calls for the insertion of a constant self-energy representative of the respective pole location where the residue is taken.We stress that such a procedure is insensitive to the analytic continuation to real frequencies as the actual self-energy is still only a constant, just not the same constant in all propagators.
In short, the presented phenomenology can be seen as the replacement ω → ω + iγ, with the additional stipulation that the relaxation rate Γ of the intermediate states is allowed to be different from relaxation rate γ Γ of the long-lived quasiparticles close to the Fermi surface.We emphasize that it is certainly possible to systematically include the very same phenomenology in the density matrix approach with length gauge, but it seems to lack the intuitive interpretation available in the diagrammatic analysis.

III. ANOMALOUS ACCELERATION
In order to intuitively understand second order response in the absence of T , one might want to fall back on a semiclassical treatment of the problem.
Disregarding lifetime effects and assuming T , the semiclassical result for the non-linear conductivity is stated as [8,26] where the relaxation rate can be reintroduced after the fact by the replacement ω → ω + iγ.The Abelian Berry curvature is defined as Ω ab mm = ∂ a r b mm − ∂ b r a mm .However, this phenomenology breaks down once the result from perturbation theory is considered.Firstly, the usual semiclassical derivation cannot account for processes which are instrinsically finite-ω, but non-negligible after the momentum space integral.Both the shift current and the injection current are of this type.Secondly, non-linear response can involve interband transitions which are easily missed in the construction of the collision integral.The situation is complicated by the fact that finite lifetimes also modify the relation between Berry connection and velocity matrix elements according to Eq. ( 19), as we discussed before.

A. Vertices and sum rules
The calculation proceeds as follows.We first note that the derivative of the Berry connection is ∂ b r c mn = R bc mn + i[r b , r c ] mn , where 2R bc mn = i m|∂ a ∂ b n + h.c., as mentioned.This matrix element can be related to the quantum metric through 2g ab mn = 2iR ab mn − l f l (r b ml r c ln + r c ml r b ln ).The Berry dipole is given by [cf.Sec.A] where all matrix elements carry two band indices which are circled through as dictated by the commutators (i.e. [a, b] mn = l (a ml b ln − b ml a ln )).We will henceforth suppress band indices whenever possible.Furthermore for interband terms it holds that i∂ c v b = −ε∂ c r b −∆ c v b +O(Γ).These relations, while not exact for finite relaxation rates, are perturbative in relaxation rate corrections and do not introduce any shifted poles.Sometimes, it is suggested to use the substitution r a mn → iv a mn /ε nm .However, as explained earlier, such is only possible in a clean system.This is noteworthy since the decomposition of w ab mn creates, among others, terms which are very similar to diagrams (I) and (II).Going forward, we will only employ the opposite route of identification via v a mn → iε mn r a mn + O(Γ), which always remains viable.This allows to immediately and transparently characterize the role of sum rules for systems with finite lifetimes and to establish a unique procedure to translate between the different gauges.
We rewrite w ab and u abc by expanding the derivatives, The decomposition of the vertices might still appear pretty verbose, but crucially it will allow for a straightforward identification of the physical processes.Independent of this, it makes manifest that relaxation rates are accounted for perturbatively in the velocity expansion of the higher order vertices.As an aside, note that only the diagonal component of u abc was analyzed since this is the only manner in which it appears here.The off-diagonal elements of u abc needed for higher orders follows straightforwardly by continuing the expansion of ∂ a ∂ b v c for off-diagonal matrix elements.Keeping the shorthand notation, the second-order response is generally While all top level commutators carry the indices mm, they do not evaluate to zero since the Fermi-Dirac factor f m is not incorporated, i.e. the expression is not a trace.Nevertheless, thanks to the Fermi factor any top level commutators actually evaluate to zero for the diagonal pieces where a denominator ε = 0 could otherwise lead to issues.This allows the replacement of velocity matrix elements by the Berry connection in outer commutators and in particular renders the iγ in the first line harmless.
In principle, the frequencies ω 1 and ω 2 in front of the integral also acquire an imaginary part according to adiabatic switching (ω → ω + iγ), but this replacement does not lead to any further complications and will henceforth be implicitly assumed for clarity of presentation.
The usual simplifications for a gapped and time-reversal symmetric system resulting in Eq. ( 20) are implemented using the sum rules Eqs.(22,23), which present the proper extensions of the classical sum rules for non-quadratic Hamiltonians in systems with finite quasiparticle lifetimes [5,7,61].T dictates that velocity matrix elements fulfill v a mn (k) = −v a mn (−k) * , while the non-Abelian Berry connection transforms like r a mn (k) = r a mn (−k) * , and R ab mn (k) = −R ab mn (−k) * .For a gapped system the dc conductivity follows from the simultaneous expansion in small ω ω 1 , ω 2 ; we will also keep track of finite lifetimes whenever necessary.The integrand of the leading order term of size ω −2 is explicitly (δ c mn = r mm − r nn ) Here, we used that We reiterate that sum rules basically constitute expansions of the vertices w and u in terms of acceleration (w) and jerk (u).Accordingly, the sum rules contain commutators with exactly the required number of odd components under momentum inversion.For example, it is Unlike in linear response, the Drude-type term ∂ a ∂ b v c is odd under momentum inversion and equally cancels.In essence, in the presence of T , higher order vertices do not contribute at leading order in ω to σ ab;c .However, from the definition of the vertices (w ab mn = m|∂ a ∂ b H|n and u abc mm = m|∂ a ∂ b ∂ c H|m ), the same conclusion can be drawn immediately from their symmetry properties under momentum inversion, thus sidestepping the need for sum rules.

B. Finite frequency terms
The last two lines in Eq. ( 24) require more care as they contain a relaxation rate in the nested denominator, where a finite γ is not zeroed out a priori.While Eq. ( 25) suggests that both terms proportional to γ −1 cancel each other, this is not necessarily the case as the limit ω → 0 does neither commute with the k-space integration nor with the limit Γ → 0. Thus, while Eq. ( 25) correctly captures the response at ω = 0, it misses not only the injection current, but also the shift current, both are by construction finite frequency responses.Their incorporation is of course straightforward by keeping all imaginary parts in Eq. (24).Denominators containing Γ have a non-zero imaginary part for ε = ω 1,2 and lead to the shift current, while denominators with γ lead to the injection current.Sipe and Shkrebtii [5] discussed the order of limits γ → 0, ω → 0. The divergence then manifests in a finite current injection, dJ/ dt = 0. Another alternative is to discuss the dc case as the limit ω → 0, followed by γ → 0 [8].In this context it is important to keep in mind that any finite quasiparticle lifetime will limit the phenomenon of current injection to irradiation times shorter than the lifetime.At longer times, the current instead saturates at a large but finite current.For simplicity, we continue to call a current proportional to γ −1 injection current, even when keeping γ finite.
The injection current is due to the conductivity In the clean limit the last term simplifies to the injection current in the form written in Eq. ( 8).In particular, for applied frequencies ω Γ, D assumes the form of δ-functions.This means that current injection will occur invariably for any combination of limits ω → 0, γ → 0, ω → 0. However, completing the expansion under the assumption that ω Γ, the result is instead This expression is finite as long as Γ/γ < ∞, meaning that the effect of current injection ceases in the strict zero frequency limit.In other words, it requires not only timereversal breaking, but also well defined quasiparticles to prevent the cancellation of the time-reversed paths.This effect is reminiscent of antilocalization in a disordered system with spin-orbit coupling.In systems with antilocalization, pairs of counterpropagating paths will not interfere destructively due to the spin-momentum locking, leading to an increased value of the linear conductivity.Upon removing spin-orbit coupling, both paths acquire opposite Aharonov-Bohm phases and the destructive interference between them leads to a decrease in conductivity.Notably, weak antilocalization persists in the strict adiabatic limit.In contrast, the effect of current injection discussed here relies on applied frequencies which exceed the inverse lifetime of excited states (ω Γ), and is thus intrinsically a dynamical effect, which does not extend to the adiabatic limit where all frequencies approach zero uniformly.Current injection could therefore be described as an effect of dynamical antilocalization.It is tied to the appearance of fermion loops with more than two vertices, meaning similar effects should be present in non-linear response functions of any order, given some suitable order of limits for the external frequencies.
Using the expression for the clean case, the δ-functions allow the substitution ω → −ε mn in the global coefficient.The result coincides with the corresponding expression in Eq. ( 8), Then the following suggestive formulation of the commutator is possible: where we used ∂ c v b = iε∂ c r b + i∆ c r b .We therefore interpret the injection current as a mismatch between v a ∂ c r b which contains the derivative of the Berry connection ∂ c r b , and r a ∂ c v b due to the kinematic acceleration ∂ c v b , each of them integrated over the manifold where ε mn = ±ω.This reinforces the dynamical nature of the injection current, which appears due to the different accelerations acting on a band wavefunction and on a quasiparticle.It is not possible to formulate either acceleration purely in terms of the Berry curvature, with the exception of a two band model, where ∆ c mn can trivially be moved outside of the commutator [19].
The same treatment can be applied to the shift current by analyzing resonances for Γ → 0, which yields This expression for the shift current is identical to the established formulation in terms of a shift vector [2].For linear polarized light it is once again helpful to insert the band energy differences (i.e.replace ω → ±ε mn ), to obtain for finite relaxation rates where we used The shift current can thus be viewed as a geometric effect with two sources, the quantum metric and the motion due to the Berry dipole, but with the intraband acceleration projected out.For an explanation of this second part, we recall that δ c r a appears as a consequence of replacing the velocity v c in Eq. ( 24) by the Berry connection, compensating the overcounting that happens with respect to the intraband motion.This intraband motion is actually diabatic and is exactly what gives rise to the injection current; consequently the projected quantity [r c , r a ]+δ c r a comprises only the remaining adiabatic processes.We emphasize that the obtained expression for the shift current is in agreement with the previous literature (e.g.[5]).In summary, borrowing from the language of antilocalization, current injection occurs because there is a mismatch between adiabatic motion and finite frequency motion, which by the order of limits which defines optical response cannot cancel.We repeat that T enforces that current injection appears only in response to circular polarized light.The shift current is in turn the sum of all remaining finite frequency processes, except for the type of intraband motion which already gives rise to the injection current.The fundamental difference between shift and injection current can also be discussed in terms of the units of the involved integrands.While the matrix elements in the shift current, Eq. ( 35) naturally form a quantity with the units of a Berry dipole, this is not true for the injection current, because there is an additional relaxation time in the denominator.This distinction will also prove useful in the semiclassical expressions.

C. Anomalous semiclassical terms
Finally, we derive the semiclassical expansion in the absence of T .To repeat, in this case Eqs.(22,23) are no longer helpful in canceling even powers in the ωdependence.Up to order ω −1 , the regular terms are As explained, injection and shift current do not appear here as they are not amenable to an expansion in small ω.Also, the leading correction for finite lifetimes is given by the replacement ω → ω + iγ.From the systematics of expression Eq. ( 24), one can immediately deduce that only the term of order ω −1 can have a purely geometric origin.All other orders necessarily depend on velocity matrix elements or the dispersion relation.Importantly, for T we recover immediately that S −2 = 0, and S −1 becomes an only function of the Berry dipole ∂ c Ω ab .The second order response, Eq. ( 24) thus decomposes into two pieces given by Eq. (29,35) with a singular frequency dependence and the regular terms, Eq. (37).Compared to previous results, the reported formulas have no further derivatives which act on the propagators, they explicate the frequency dependence and they phrase the response in terms of four physical processes, each of which is, in principle, measurable independently.As an aside, we point out that the most often used formula for shift and injection current, Eq. ( 8), does not lend itself to a straightforward generalization for finite lifetimes, as the replacement of applied frequencies by band energy differences is only possible in the limit ω γ.Instead, the correct starting point in the length gauge are expressions containing the covariant derivative, which can be found for example in [7,47].
In the absence of a gap, the semiclassical expansion, Eq. ( 37) for small ω becomes questionable.This is, however, not a problem for the wide range of Dirac and Weyl semimetals which have undertilted cones: The regions in k-space in the direct vicinity of the Dirac points do not contribute significantly to the non-linear conductivity, rendering them irrelevant.In contrast, in type II Weyl semimetals this does not hold, hence they exhibit a very different phenomenology [64].

D. Conductivity tensor
At this point it is helpful to recall that sometimes several components of the non-linear conductivity σ ab;c are combined for a given current response, depending on the applied electric fields.For example, the BPVE for linear polarized light with polarization direction along the x, y-diagonal is given by In the same fashion, light with anticlockwise polarization in the x, y-plane leads to the current We will employ some of these properties in the next section to recast the discussion in a more familiar fashion in terms of the symmetries of real and imaginary parts of the matrix elements.

IV. DC CURRENT RESPONSE
In the following we explore the dc-current creation from a single Weyl cone and relate it to the general phenomenology we explained in the previous section.Of particular interest will of course be the injection current, but we also make some comments regarding the remaining contributions which combine to the shift current.
Since the Abelian Berry curvature is strictly zero in a space-time inversion (combined symmetry by the inversion symmetry and T , PT ) symmetric material, the contribution from the Berry dipole vanishes.We utilize this special case to investigate the dc-response when it is dominated by acceleration terms unrelated to the Berry curvature.
To repeat, the second-order optical response is given by the six diagrams depicted in Fig. 1.The dc-current response is σ ab;c (0, ω, −ω) ≡ σ ab;c dc .Diagrams (III)-(VI) contain the contribution from curvature effects in the intraband motion, they do not contribute to the injection current and will therefore be disregarded in this section.Additionally, for the benefit of familiarity, in this section we will fall back on the standard notation in band indices instead of the very compact notation of the previous section.Without loss of generality, consider in the following a coordinate system which is aligned with the direction of the incident linear polarized light.Then for the circular polarization, the indices a and b of the electric field are dissimilar while for the linear polarized light it is a = b.For any other choice of coordinates, the following statements still hold for the respective current response, but not individually for the matrix elements.It follows in either case that Furthermore, where Note that the intraband relaxation rate γ is usually significantly smaller than the interband relaxation rate Γ.The formulation as real and imaginary part follows form the observation that the incoming electric fields either fulfill polarization).We will occasionally use the shorthand Σ ab;c ± (ω) with σ abc lin = Re Σ ab;c + (ω) and σ abc circ = Im Σ ab;c − (ω) A. Cancellations from time reversal symmetry In the presence of T , energy eigenvalues are inversion symmetric in the Brillouin zone, ε n (k) = ε n (−k) and velocity matrix elements fulfill v a nm (k) = −v a nm * (−k).In this case, both injection current and shift current benefit from cancellations between k and −k.The important term to keep track off is ε mn −iγ, which in the limit γ → 0 is manifestly divergent for m = n.We therefore proceed to write out the real and imaginary parts of the energy denominators D −1 mn = ε mn − iγ and D −1 nl (ω) = Ω − ε nl − iΓ.Since the nonlinear dc-response is proportional to E a (−ω)E b (ω) one can replace a ↔ b and ω ↔ −ω in the expression for A, meaning that for an inversion symmetric Brillouin zone the second order response can be written as proportional to where s = +1 (s = −1) for linear (circular) polarized light.For the various combination possibilities of real and imaginary parts of J lin , this implies In Eqs.(46,49) we made use of the fact that a = b for linear polarized light, which implies Re V 1 = 0. Eq. ( 47) follows from the fact that Im D mn = 0 requires m = n, which entails Im V 1 = 0, and Eq. ( 51) is zero due to the antisymmetry in the band summation for m = n.

B. Space-time inversion symmetry
In the absence of T , the aforementioned cancellations do no longer hold and a comparable behavior (injection current) is expected for both linear and circular polarization.Using the conventional argumentation, without T any combination of real and imaginary parts of the denominators is potentially nonzero and has to be reconsidered.This is of course nonwithstanding crystal symmetries, which can still render certain components zero.As we discussed, the well-known formula for the shift current and linear polarization, is no longer the most relevant part of the response function once T is broken.This term originates from Eq. ( 48) only and evaluates to zero for m = n.Instead, resonances involving m = n provide the dominant contribution to the BPVE.It is then useful to separate two-band contributions containing Im D mn = γ −1 δ mn from three-band ones with Re D mn = ε −1 mn (1 − δ nm ), the latter of which are suppressed by the energy denominator ε mn .
In the presence of PT -symmetry, only the real part of the product of momentum matrix elements N abc nml contributes.Imposing PT and focusing on the dominant two-band contributions, the response therefore features the following term Note that for a T material, only the imaginary part of the numerator N nml contributes to the dc-current, while in the presence of PT -symmetry only the real part does.In this sense, both are complementary.Since T is broken, a symmetrization with respect to k is no longer useful.Instead, one should consider each Dirac cone separately.
It was previously noted that semimetals have zero shift current if the cone has a linear dispersion [18,64].This is related to the triple product of velocities, which is antisymmetric on a radial shell which concentrically encloses a perfect Dirac point.
To be concrete, suppose a tight binding model with PTsymmetry, constituted of two pairs of degenerate bands.Given a trivial implementation of inversion symmetry, the most general PT -symmetric four-band Hamiltonian is [65] where σ and τ are Pauli matrizes of spin and orbital degrees.If the orbital make-up requires a nontrivial implementation of inversion, some spatial indices will be permuted.Due to the high symmetry in this four-band model, for g x0 = g yx = g yy = 0, the basis can be chosen pairwise such that one can express dispersion and momentum matrix elements referring to only one upper and one lower band index.In the general case one obtains not one but three structurally identical pieces corresponding to the three different interband transitions.Inserting the degeneracies into the general expression, and using that ω nn = 0, ω mn = −ω nm , v a nm = v a mn * , one obtains for the non-linear dc conductivity explicitly where the upper (lower) band is labeled u (d).The dots indicate that there are analogous terms including the other three pairings of bands (involving the T -partners ū, d), where The problem thus decomposes in four pieces which somewhat resembles the two-band case when T is preserved [27].However, now there are both two-band and three-band terms.

C. Single cone
Concentrating on one band touching point at a time, the local two-band Hamiltonian corresponding to band (d, u) is H(k) = g i σ i , with a generic momentum dependencies g i = v F (k i + α ijl k j k l ); analogous expressions with dissimilar coefficients exist for the other combinations of bands at the remaining Weyl points.Keeping only the two-band expression for linear polarization, the injection current becomes The velocity matrix elements are given by Re Expanding in the non-linearity α, the zeroth term vanishes, as expected.At the next order, since i = j in Eq (59), the only contribution is (Γ → 0) so the two cases of linear polarization are (61) For finite interband relaxation rates (Γ = 0), there is a crossover at frequency ω c ∼ (ΓΛ 2 ) 1/3 , with Λ the bandwidth of the system.Below this ω c , σ aa;c dc acquires a residual term of the size ∼ ΓΛ 2 /γω 2 .
Barring any symmetry cancellations, for transitions near the Dirac point we thus obtain the estimate where c 1 , c 2 , c 3 are coefficients depending on the nonlinearity in the Dirac cone dispersion.The non-identical Dirac cone close to −k will contribute a similar term with numerically different coefficients, but this does not change the the qualitative behavior.Additionally, the values of c 1 and c 3 are not exclusively determined by diagrams (I) and (II) but will also receive a contribution from the remaining second-order diagrams.In summary, in case that both T and inversion are preserved, the second order dc-response will vanish.With T present but P being absent, what matters is the imaginary part of the product of three velocities (N abc nml ), with P present but without T , both the real and imaginary part of N abc nml matters, and in the absence of both P and T but with T -symmetry only the real part of N abc nml remains.For the T -symmetric case, weak T breaking will therefore result in a strong suppression of the non-linear dc-conductivity due to the nearly complete cancellation of the contributions from k with the ones from −k.A strong response is therefore expected from intrinsically T -breaking materials where the T is not broken perturbatively.At leading order the size of the injection current is dependent on the band structure near the Dirac cone through the combination α/γ, which means that a Dirac cone with a large non-linear components of type α cac or α caa appears to be the most promising candidate for measuring a large dc-current.
Regarding the contributions from diagrams (III)-(VI) which we safely discarded, the situation becomes quite different if T is preserved, rendering diagrams (I) and (II) of size O(1) instead of O(γ −1 ).Then, all diagrams have a similar sized contribution to the dc-current [27].

D. Minimal four band model
We corroborate the findings of the previous section with an explicit calculation of the BPVE in a T -symmetric four band model.As a comparison, a similar four-band model exhibiting T with comparable dispersion and density of states is chosen.While both models are somewhat synthetic, they are a suitable representative of the BPVE which arises from two-photon processes in a Dirac cone close to the nodal point.Of course, at higher frequency resonances will appear which are particular to the chosen band structure.
As tight binding Hamiltonian we choose, The parameters m y , m z produce anisotropy, ∆ breaks inversion symmetry.δ induces time-reversal and inversion breaking, such that while the Hamiltonian does not commute with either P or T , their combined operation PT is a symmetry.The spectrum is then of two doubly-degenerate bands.In the following we compare the non-linear conductivity between the choice of parameters ∆ = 0, δ = 0, For improved visibility, we show the product σ yy;z ω 2 .For PT , the non-linear dc-conductivity approaches a ω −2 -dependence for small frequencies which is independent of the quasiparticle lifetime γ.At frequencies above a threshold proportional to Γ, the current creation is only capped by the lifetime 1/γ.For T , the current vanishes linearly with ω 1 for small frequencies, save for a residual term of size γω −2 for the smallest frequencies.At frequencies comparable to the bandwidth, the spectrum is sensitive to nesting effects, and thus to the band structure.Note that the apparent erratic behavior for these large values of ω is due to the logarithmic scale, not numerical artefacts.
Parameter values are the same as in Fig. 3.
implementing T and the choice ∆ = 0, δ = 0 for PTsymmetry.The band structures have the form shown in Fig. 3.Both cases lead to a very similar total density of states (DOS) (Fig. 3c), with a clear separation between the low energy Dirac cone physics and the high energy Lifshitz point where the cones connect.The Hamiltonian, Eq. ( 65) is generic in the sense that it encodes the minimal number of four Weyl cones in case of T and the minimal number of two doubly degenerate Dirac cones in the PT case.The momentum dependencies are not generic, they are chosen in reference to a minimal model which was proposed for the time-reversal symmetric material TaAs, which has a cubic Brillouin zone.In Sec.B, the steps are listed which lead to Eq. ( 65), suffice to say here that the only nontrivial part is to keep enough structure in the momentum dependencies such that there are no unwanted additional symmetries like mirrors present.There is one specialty added to the PTsymmetric dispersion which is to not include a term which would correspond to a non-linearity in the dispersion of the type α zyy in the vicinity of the Dirac points.This allows us to verify the claim of Eq. ( 63) that the injection current is sourced from exactly one non-linear dispersion term at next-to-leading order.In other words, the expected injection current for our PT -symmetric example for the dc-conductivity σ yy;z has size O(ω 2 /γ) instead of O(ω/γ).For the numerical evaluation we repeat that diagrams (III)-(VI) are negligibly small, so that it is enough to consider diagrams (I) and (II), with varying but finite relaxation rates γ = Γ.The second order dc-conductivity σ yy;z is shown in Fig. 4 as a function of frequency.As expected, the injection current persists for large frequencies, crossing over into a shift current for small ω.We emphasize the difference in magnitude of roughly 10 2 between T and PT , in spite of the very comparable band structure, bandwidth, and relaxation rates used for both models.

V. CONCLUSIONS
For many years, studies of the bulk photovoltaic effect were restricted to the evaluation of the complicated perturbative expressions.Here, we used the diagrammatic approach to describe second-order optical response for materials without time-reversal symmetry and also with finite relaxation times in terms of four distinct physical processes.The central element was the analysis of the anomalous acceleration, which appears naturally in the velocity gauge.The resulting analytical structure (e.g.Eqs. 24, 29 and 35) of the expressions is substantially more amenable to analytical manipulations, but can also be directly implemented in numerics.
We derived explicit sum rules for non-quadratic Hamiltonians and demonstrated that their cancellation properties are equivalently accounted for in the diagrammatic formulation of the velocity gauge.
Contrary to the semiclassical intuition, the presence of a Berry dipole is not a necessary requirement for a large BPVE.Instead, the topological content of the secondorder response is encoded in expressions which involve both the antisymmetric and the symmetric part of the quantum geometric tensor.In particular, we showed that the shift current can be written as the sum of a piece containing the quantum metric and certain elements of the Berry dipole.The shift current is therefore of topological origin, but not in the usual sense that it would be proportional to the Berry curvature.Also, it does not evaluate to a quantized value in a topological insulator.In agreement with earlier work, we recover the Berry dipole in the semiclassical expansion, but alongside find a number of other terms which only vanish in the presence of T .This means that the semiclassical current response is of topological origin only in time-reversal symmetric compounds.From the same analysis, we conclude that the injection current is generically not of topological origin, instead being attributed to dynamical antilocalization between counter-propagating momentum space trajectories.
We believe that the same analysis is feasible for third order conductivities, with a larger number of physically distinct processes emerging for the matrix elements involving third derivatives ('jerk').
The principles outlined here should prove helpful in the search for materials with tailor-made non-linear optical response.For example, we outlined that current injection is enhanced in semimetals with broken time-reversal symmetry when they exhibit a large non-linear term in their dispersion near the Dirac cone.

Figure 1 .
Figure 1.Second order contributions to the optical dc-response.Straight lines with arrows denote quasiparticles, a straight plain line connects to a current operator and the wavy lines are couplings to the electromagnetic field.Note that diagrams (I) and (II) and (III) and (IV) are variants of each other.

Figure 2 .
Figure2.Phenomenology of on-shell processes.the diagram on the left for dc response contains two propagators at the same frequency, while the second harmonic generation involves propagators at three different frequencies.In the language of on-shell/off-shell processes, only for the diagram on the left one encounters a pole which is restricted to intraband relaxation rates near the Fermi surface.

Figure 4 .
Figure 4. Comparison of the modulus of σ yy;z (0, ω, −ω) between the PT -symmetric case (a) and the one with T (b).For improved visibility, we show the product σ yy;z ω 2 .For PT , the non-linear dc-conductivity approaches a ω −2 -dependence for small frequencies which is independent of the quasiparticle lifetime γ.At frequencies above a threshold proportional to Γ, the current creation is only capped by the lifetime 1/γ.For T , the current vanishes linearly with ω 1 for small frequencies, save for a residual term of size γω −2 for the smallest frequencies.At frequencies comparable to the bandwidth, the spectrum is sensitive to nesting effects, and thus to the band structure.Note that the apparent erratic behavior for these large values of ω is due to the logarithmic scale, not numerical artefacts.Parameter values are the same as in Fig.3.