Gauge invariance and Ward identities in nonlinear response theory
Habib Rostami, Mikhail I. Katsnelson, Giovanni Vignale, Marco Polini
GGauge invariance and Ward identities in nonlinearresponse theory
Habib Rostami a , Mikhail I. Katsnelson b , Giovanni Vignale c , Marco Polini d,e,f a Nordita, KTH Royal Institute of Technology and Stockholm University, Stockholm SE-10691, Sweden b Institute for Molecules and Materials, Radboud University, Heyendaalseweg 135, 6525 AJ,Nijmegen, The Netherlands c Department of Physics and Astronomy, University of Missouri, Columbia, Missouri65211, USA d Dipartimento di Fisica dell’Universit`a di Pisa, Largo Bruno Pontecorvo 3, I-56127Pisa, Italy e Istituto Italiano di Tecnologia, Graphene Labs, Via Morego 30, I-16163 Genova, Italy f School of Physics & Astronomy, University of Manchester, Oxford Road, Manchester M139PL, United Kingdom
Abstract
We present a formal analysis of nonlinear response functions in terms of correla-tion functions in real- and imaginary-time domains. In particular, we show thatcausal nonlinear response functions, expressed in terms of nested commutatorsin real time, can be obtained from the analytic continuation of time-orderedresponse functions, which are more easily amenable to diagrammatic calcula-tion. This generalizes the well-known result of linear response theory. We thenuse gauge invariance arguments to derive exact relations between second-orderresponse functions in density and current channels. These identities, which arenon-perturbative in the strength of inter-particle interactions, allow us to estab-lish exact connections between nonlinear optics calculations done in differentelectromagnetic gauges.
Keywords:
Nonlinear optics; Gauge invariance; Ward identities
1. Introduction
Nonlinear optical phenomena [1–4], arise from the response of an electronicsystem to intense electromagnetic fields. Going beyond linear response, manyinteresting optical effects can occur including harmonic generation, wave mix-ing, and saturable absorption [3]. As in the case of linear response theory,calculations of nonlinear response functions to an electromagnetic field can becarried out in different electromagnetic gauges such as the scalar potential gauge(in which the scalar potential Φ is finite, while the vector potential is A = 0) Email address: [email protected] (Habib Rostami)
Preprint submitted to Elsevier February 9, 2021 a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b nd the vector potential gauge (Φ = 0 , A (cid:54) = ). Nonlinear response functionsare usually expressed in terms of complicated and lengthy relations, which arecumbersome to evaluate [2–4]. Therefore, one often needs to introduce someapproximations in order to handle technical difficulties. For example, in theso-called “optical” limit, the photon momentum is neglected and one needs tostudy only the local response to a homogenous time-dependent electric field.This treatment is often called electric-dipole approximation [2, 3]. In this ap-proximation, the field-particle interaction can be written in two major gauges: H Φ = H ( p ) − e E ( t ) · r or H A = H ( p + e A ( t )), where H is the Hamiltonianin the absence of radiation, and r and p ≡ − i (cid:126) ∇ r stand for the position andmomentum operators, respectively. The equivalence of these two gauges for anarbitrary physical observable in the electric-dipole approximation is well dis-cussed in the literature [5–9].For optical properties of insulators, we mainly need to consider inter-bandtransitions while the intra-band contribution is more relevant in metals andsemimetals. To calculate the inter-band optical response, it is convenient toutilize the vector-potential gauge [10] in which one can simply drop the photonmomentum from the very beginning of the calculation. Apart from pure intra-and inter-band transitions, an extra class of transitions emerges only in theanalysis of nonlinear response functions containing n -point correlation functionswith n ≥
3. This new contribution originates from mixed transitions whichcontain both intra- and inter-band processes [11, 12]. Employing a genericdensity gauge, Φ( r , t ), is more suitable for evaluating the contributions due toall intra-band, inter-band and mixed transitions in the optical limit [11–13].This is because only a single non-interacting Feynman diagram is required tobe evaluated in the density gauge at any order of perturbation theory, while inthe vector potential gauge there are two, four, and eight diagrams that needto be evaluated for the linear, second-order and third-order response functions,respectively [14].In nonlinear response theory, the gauge choice is crucial when going beyondthe electric-dipole approximation. If one keeps the photon momentum up to lin-ear order in the nonlinear conductivity, it is possible to capture nonlocal effectsin the electric-quadrupole and magnetic-dipole approximations [13, 15]. Thesenew contributions are very important for the second-order nonlinear response ofinversion symmetric systems [3, 16, 17]. Note that in the scalar potential gauge,one loses the transverse contribution to the current, which originates from themagnetic-dipole moment [15]. In fact, as discussed in this Article, in order tocalculate the nonlocal current it is preferable to perform the calculation in thevector potential gauge by considering a spatially inhomogeneous vector field, A ( r , t ).As just reviewed, most of the available literature on the issue of gauge invari-ance in nonlinear response theory is devoted to exploring the gauge choice (thetwo specific aforementioned gauges in the electric-dipole approximation) in non-interacting atoms [18] and crystalline solids [8]. To the best of our knowledge, ageneric analysis of gauge invariance for a spatially inhomogeneous external fieldin nonlinear response theory as applied to interacting electron systems is still2nexplored and it is one of the main motivations of this work. Our analysis ofnonlinear response functions covers all local (electric-dipole) and nonlocal (otherelectric and magnetic-multipoles) effects and can be seen as a generalization ofthe gauge invariance analysis in the context of linear response theory, which istextbook material [19–21]. Nonlocal response functions are also of interest whenthe momentum exchanged between photons (external field) and charged parti-cles cannot be ignored. This is the case of the photon-drag effect [22–24] andof other effects where the photon momentum needs to be taken into account.For example, our analysis may be relevant to analyze experimental results ofspectroscopy based on free-electron lasers [25, 26].Gauge invariance, i.e. the independence of physical observables on the elec-tromagnetic gauge one chooses for calculations, imposes severe constraints onthe theory of linear response [19–21]. Indeed, 1) the fact that a static and purelylongitudinal vector potential cannot produce any physical current implies thevanishing of the longitudinal current response at any finite wave vector q ; 2)similarly, the fact that a static and quasi-homogeneous vector potential can-not produce any physical current implies the vanishing of the longitudinal and transverse current-current response functions for q →
0. More generally, gaugeinvariance imposes precise relationships between linear response functions toscalar and vector potentials.In this Article, we first present a formal analysis of the nonlinear responsetheory in the real and imaginary time domains. We show that causal nonlinearresponse functions, expressed in terms of nested commutators in real time, canbe obtained from the analytic continuation of time-ordered response functions,which are more easily amenable to diagrammatic calculation. This generalizesthe well-known result of linear response theory.We then provide a theoretical study of the impact of gauge invariance onthe second-order response functions. We finally report Ward identities whichmust be fulfilled in order to ensure gauge invariance at any order of perturba-tion theory. In Section 2, we report an explicit definition of nonlinear responsefunctions. In Section 3, we discuss the spectral representation of second-order re-sponse functions in real- and imaginary-time domains. In Section 4, we considera specific kind of external field, i.e. electromagnetic radiation treated classically.In Sections 5 and 6, we explicitly discuss gauge invariance for linear and second-order response functions, respectively. Finally, in Section 7, we present a set ofWard identities for nonlinear response functions.
2. Nonlinear response theory in real time
We consider a many-body system in thermal equilibrium which is describedby a Hamiltonian denoted by ˆ H . We then turn on an external field, which canbe modelled by a field-particle interaction term denoted by ˆ V . We use the ˆ S -matrix approach [27] to consider the effect of the interaction in a perturbativeway. In this approach, one has an adiabatic time evolution of the wave functionfrom an unperturbed state at time t = −∞ , i.e. | ψ (cid:105) , to the perturbed state at3ime t . Therefore, we have | ψ ( t ) (cid:105) = ˆ S ( t, −∞ ) | ψ (cid:105) where the ˆ S -matrix is givenby (setting (cid:126) = 1 for the sake of simplicity)ˆ S ( t, −∞ ) = T exp (cid:18) − i (cid:90) t −∞ ˆ V ( t (cid:48) ) dt (cid:48) (cid:19) . (1)Here, T stands for the time-ordering operation and ˆ V ( t ) is the perturbative partof the Hamiltonian in the interaction representation, i.e. ˆ V ( t ) = e it ˆ H ˆ V e − it ˆ H . Ina similar way, any observable evolves in the interaction picture asˆ A ( t ) = e it ˆ H ˆ Ae − it ˆ H . (2)Here, ˆ A denotes the second-quantized representation of the A operator. Theexpectation value of ˆ A at time t is given by A ( t ) = (cid:104) ψ ( t ) | ˆ A ( t ) | ψ ( t ) (cid:105) = (cid:104) ψ | ˆ S † ( t, −∞ ) ˆ A ( t ) ˆ S ( t, −∞ ) | ψ (cid:105) . (3)Thanks to the Dyson expansion, we have (see Appendix A) A ( t ) = (cid:68) ˆ A ( t ) (cid:69) + i (cid:90) ∞−∞ dτ Θ( τ ) (cid:68)(cid:104) ˆ V ( t − τ ) , ˆ A ( t ) (cid:105)(cid:69) + i (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ Θ( τ − τ )Θ( τ ) (cid:68)(cid:104) ˆ V ( t − τ ) , (cid:104) ˆ V ( t − τ ) , ˆ A ( t ) (cid:105)(cid:105)(cid:69) + i (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ Θ( τ − τ )Θ( τ − τ )Θ( τ ) × (cid:68)(cid:104) ˆ V ( t − τ ) , (cid:104) ˆ V ( t − τ ) , (cid:104) ˆ V ( t − τ ) , ˆ A ( t ) (cid:105)(cid:105)(cid:105)(cid:69) + . . . , (4)where Θ( τ ) is the Heaviside function. The field-particle interaction can beformally written as followsˆ V ( t ) = (cid:88) i ˆ B i ( t ) F i ( t ) + (cid:88) ij ˆ C ij ( t ) F i ( t ) F j ( t ) + (cid:88) ijk ˆ D ijk ( t ) F i ( t ) F j ( t ) F k ( t ) + . . . , (5)where ˆ B i ( t ) indicates linear coupling while ˆ C ij ( t ) and ˆ D ijk ( t ) stand for the non-linear couplings. We can visualize these linear and nonlinear coupling terms asvertices in the diagrams shown in Fig. 1. For the case of light-matter inter-actions, we term the B , C , and D vertices as single-photon, two-photon, andthree-phonon vertices, respectively. Here, we first assume a linear coupling viaˆ B i to the external fields F i and later generalize our formal theory to the non-linear couplings ˆ C ij ( t ) and ˆ D ijk ( t ). Note that vertex couplings like ˆ C ij ( t ) andˆ D ijk ( t ) only emerge in the non-relativistic limit of a low-energy theory becauseof the nonlinear energy dispersion of the band Hamiltonian ( (cid:15) k ∼ (cid:80) n ≥ (cid:15) n k n )and they do not exist in quantum electrodynamics [28]. After neglecting non-4 B i ˆ C ij ˆ D ijk Figure 1: Diagrammatic representation of the linear and nonlinear couplings as vertices (soliddots) coupled to external fields (dashed lines). In the case of electromagnetic perturbations,we call the B , C , and D vertices as single-photon, two-photon and three-phonon vertices,respectively. linear coupling terms and plugging Eq. (5) in Eq. (4), we obtain A ( t ) = (cid:68) ˆ A ( t ) (cid:69) + (cid:88) i (cid:90) ∞−∞ dτ χ AB i ( t ; τ ) F i ( t − τ )+ (cid:88) ij (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ χ AB i B j ( t ; τ , τ ) F i ( t − τ ) F j ( t − τ )+ (cid:88) ijk (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ χ AB i B j B k ( t ; τ , τ , τ ) × F i ( t − τ ) F j ( t − τ ) F k ( t − τ ) + . . . . (6)The previous expression leads us to define the linear retarded response func-tions [21] χ AB i ( t ; τ ) = i Θ( τ ) (cid:68)(cid:104) ˆ B i ( t − τ ) , ˆ A ( t ) (cid:105)(cid:69) , (7)the second-order retarded response functions χ AB i B j ( t ; τ , τ ) = i (cid:48) (cid:88) P Θ( τ − τ )Θ( τ ) (cid:68)(cid:104) ˆ B j ( t − τ ) , (cid:104) ˆ B i ( t − τ ) , ˆ A ( t ) (cid:105)(cid:105)(cid:69) , (8)and the third-order ones χ AB i B j B k ( t ; τ , τ , τ ) = i (cid:48) (cid:88) P Θ( τ − τ )Θ( τ − τ )Θ( τ ) × (cid:68)(cid:104) ˆ B k ( t − τ ) , (cid:104) ˆ B j ( t − τ ) , (cid:104) ˆ B i ( t − τ ) , ˆ A ( t ) (cid:105)(cid:105)(cid:105)(cid:69) . (9)Of course, one can simply extend these definitions to higher-order responsefunctions. The symbol (cid:80) (cid:48)P in Eqs. (8) and (9) is there to ensure a permutationsymmetry. Since the quantity Π mj =1 F j ( t − τ j ) in Eq. (6) is symmetric with respectto the permutation between each pair of dummy variables, i.e. i ≡ ( F i , τ i ) and j ≡ ( F j , τ j ), we expect a permutation symmetry for the m -th order responsefunction. For instance, following Ref. [3], one can decompose the m -th orderresponse function in the sum of symmetric and anti-symmetric contributions,5.e. χ = χ s + χ a with χ s AB B ...B m ( t ; τ , τ , . . . , τ m ) = 12 (cid:110) χ AB B ...B m ( t ; τ , τ , . . . , τ m )+ χ AB B ...B m ( t ; τ , τ , . . . , τ m ) (cid:111) (10)and χ a AB B ...B m ( t ; τ , τ , . . . , τ m ) = 12 (cid:110) χ AB B ...B m ( t ; τ , τ , . . . , τ m ) − χ AB B ...B m ( t ; τ , τ , . . . , τ m ) (cid:111) . (11)It is obvious that χ s ( χ a ) is symmetric (anti-symmetric) with respect to thepermutation between ( B , τ ) and ( B , τ ). Since the term Π mj =1 F j ( t − τ j ) issymmetric for all orders of permutation, the anti-symmetric part plays no rolein determining the expectation value of ˆ A . The only non-vanishing contributionto the latter originates from completely symmetric response functions. This isthe reason why in the definition of the m -th order nonlinear response functionwe have introduced the sum (cid:80) (cid:48)P = m ! (cid:80) P , where (cid:80) P stands for the sum overall permutations among the dummy variables ( B i , τ i ). This kind of symmetryis usually called intrinsic permutation symmetry [3].By using the cyclic properties of the trace, we can then eliminate “ t ” inEqs. (7), (8), and (9): χ AB i ( τ ) = i Θ( τ ) (cid:68)(cid:104) ˆ B i ( − τ ) , ˆ A (0) (cid:105)(cid:69) , (12) χ AB i B j ( τ , τ ) = i (cid:48) (cid:88) P Θ( τ − τ )Θ( τ ) (cid:68)(cid:104) ˆ B j ( − τ ) , (cid:104) ˆ B i ( − τ ) , ˆ A (0) (cid:105)(cid:105)(cid:69) , (13) χ AB i B j B k ( τ , τ , τ ) = i (cid:48) (cid:88) P Θ( τ − τ )Θ( τ − τ )Θ( τ ) × (cid:68)(cid:104) ˆ B k ( − τ ) , (cid:104) ˆ B j ( − τ ) , (cid:104) ˆ B i ( − τ ) , ˆ A (0) (cid:105)(cid:105)(cid:105)(cid:69) . (14)The invariance under time translations of the previous response functions isevident. They just depend on the time differences τ i = t − t i . This can bepowerfully used by Fourier transforming to the frequency domain: χ AB i ( ω ) = (cid:90) ∞−∞ dτ χ AB i ( τ ) e i ( ω + iη ) τ , (15) χ AB i B j ( ω , ω ) = (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ χ AB i B j ( τ , τ ) e i ( ω Σ + i (cid:80) i η i ) τ i , (16) χ AB i B j B k ( ω , ω , ω ) = (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ χ AB i B j B k ( τ , τ , τ ) × e i ( ω Σ + i (cid:80) i η i ) τ i . (17)6 B m ˆ B ˆ B ˆ B ˆ B ˆ B ˆ A. . . − τ − τ − τ − τ − τ − τ m (a) ˆ B m ˆ B ˆ B ˆ B ˆ B ˆ B ˆ A. . . ω Σ − ω − ω − ω − ω − ω − ω m (b) Figure 2: Feynman diagram for the m -th order response function in the time, panel (a),and frequency, panel (b), domains. These diagrams just stand for one of the m ! possiblepermutations among the dummy variables. Notice that ω Σ ≡ (cid:80) mi =1 ω i . Here, ω Σ = (cid:80) i ω i . As usual [21], we have assumed that all frequencies containan infinitesimal positive imaginary part η i → + . This stems from the needto fit periodic perturbations F i ( t ) = F i ( ω ) e − iωt e − ηt + c . c . into the responsetheory formalism, making sure that these vanish in the remote past. In thefrequency domain, the intrinsic permutation symmetry is among the ( B i , ω i )pairs. The thermal average in Eqs. (12), (13), and (14) is defined as usual [21],i.e. (cid:104) ˆ O(cid:105) = Tr[ˆ ρ ˆ O ] = (cid:80) λ P λ (cid:104) λ | ˆ O| λ (cid:105) where ˆ ρ is the density matrix, P λ = (cid:104) λ | ˆ ρ | λ (cid:105) = e − βE λ / Z , Z = (cid:80) λ e − βE λ is the partition function, and β = 1 / ( k B T )with T the temperature. Notice that | λ (cid:105) and E λ are the exact eigenstates andeigenvalues of the many-body Hamiltonian ˆ H in the absence of radiation. Thesame assumptions that are made on P λ in the linear-response Kubo formalism(see footnote 11 in Chapter 3 of Ref. [21]) are assumed to hold true in thisnonlinear case too.With the definitions given in Eqs. (12), (13), and (14), we can take advantageof diagrammatic techniques to evaluate nonlinear response functions. In Fig. 2,we give a diagrammatic representation of the m -th order response function inboth time and frequency domains. The generalization of the linear-couplingtheory to the nonlinear couplings is straightforward as it is simply based on theinclusion of higher-order photon vertices in the Feynman diagrams.
3. Spectral representation and analytic continuation
In an interacting many-body system, it is not easy to handle all the diagramsfor linear and nonlinear response functions. It is often convenient to do the cal-culations with time-ordered response functions in imaginary time and then comeback to real time by an analytic continuation [27]. The imaginary-time repre-sentation of response functions, also known as “Matsubara representation”, isalso extremely useful for the finite temperature analysis. In order to derivethe connection between causal (real-time) response function and time-ordered7esponse functions in imaginary time, we proceed in three steps. First we de-rive the so-called spectral representation (also known as “Lehmann representa-tion” or “exact eigenstates representation”) for second-order response functions.Next, making use of this representation, it is shown that the second-order causalresponse functions in real time are connected to time-ordered second-order re-sponse functions in imaginary time by the well-known analytic continuationprocedure [27]. Finally, reasoning by induction, we prove that the analyticcontinuation procedure works for response functions of arbitrary order.
According to Eq. (13), the second-order retarded response function in thereal-time representation reads as follows: χ AB B ( τ , τ ) = i (cid:48) (cid:88) P Θ( τ − τ )Θ( τ ) (cid:68)(cid:104) ˆ B ( − τ ) , (cid:104) ˆ B ( − τ ) , ˆ A (0) (cid:105)(cid:105)(cid:69) = i (cid:48) (cid:88) P Θ( τ − τ )Θ( τ ) I AB B ( τ , τ ) , (18)where I AB B ( τ , τ ) = (cid:68) ˆ B ( − τ ) ˆ B ( − τ ) ˆ A (0) (cid:69) − (cid:68) ˆ B ( − τ ) ˆ A (0) ˆ B ( − τ ) (cid:69) − (cid:68) ˆ B ( − τ ) ˆ A (0) ˆ B ( − τ ) (cid:69) + (cid:68) ˆ A (0) ˆ B ( − τ ) ˆ B ( − τ ) (cid:69) . (19)Performing the thermodynamic average, we obtain I AB B ( τ , τ ) = (cid:88) λ λ λ P λ (cid:40) (cid:104) λ | ˆ B ( − τ ) | λ (cid:105)(cid:104) λ | ˆ B ( − τ ) | λ (cid:105)(cid:104) λ | ˆ A (0) | λ (cid:105)− (cid:104) λ | ˆ B ( − τ ) | λ (cid:105)(cid:104) λ | ˆ A (0) | λ (cid:105)(cid:104) λ | ˆ B ( − τ ) | λ (cid:105)− (cid:104) λ | ˆ B ( − τ ) | λ (cid:105)(cid:104) λ | ˆ A (0) | λ (cid:105)(cid:104) λ | ˆ B ( − τ ) | λ (cid:105) + (cid:104) λ | ˆ A (0) | λ (cid:105)(cid:104) λ | ˆ B ( − τ ) | λ (cid:105)(cid:104) λ | ˆ B ( − τ ) | λ (cid:105) (cid:41) . (20)We exchange λ , λ and λ dummy labels in the second term of the aboverelation in order to the collect it with the first term. We use a similar trick andcollect the last two terms together. Consequently, we arrive at I AB B ( τ , τ ) = (cid:88) λ λ λ A λ λ (cid:40) P λ λ (cid:104) λ | ˆ B ( − τ ) | λ (cid:105)(cid:104) λ | ˆ B ( − τ ) | λ (cid:105)− P λ λ (cid:104) λ | ˆ B ( − τ ) | λ (cid:105)(cid:104) λ | ˆ B ( − τ ) | λ (cid:105) (cid:41) . (21)8otice that (cid:104) λ | ˆ A (0) | λ (cid:105) = A λ λ and P λ λ = P λ − P λ . We extract thetime-dependent part as follows I AB B ( τ , τ ) = (cid:88) λ λ λ A λ λ (cid:40) P λ λ B ,λ λ B ,λ λ e iτ E λ λ e iτ E λ λ − P λ λ B ,λ λ B ,λ λ e iτ E λ λ e iτ E λ λ (cid:41) , (22)where E λ λ = E λ − E λ . Using Eq. (16), we obtain the frequency-domainrepresentation of the second-order response function: χ AB B ( ω , ω ) = i (cid:48) (cid:88) P (cid:90) ∞ dτ (cid:90) ∞ τ dτ I AB B ( τ , τ ) e i (cid:80) i ( ω i + iη i ) τ i . (23)Considering Eq. (22), we can easily perform the integrations in Eq. (23) and wereach χ AB B ( ω , ω ) = (cid:48) (cid:88) P (cid:88) λ λ λ A λ λ ω + ω + E λ λ + i [ η + η ] × (cid:40) B ,λ λ B ,λ λ P λ λ ω + E λ λ + iη − B ,λ λ B ,λ λ P λ λ ω + E λ λ + iη (cid:41) . (24)Eventually, by considering the intrinsic permutation, we arrive at the followingfinal expression for the retarded second-order response function: χ AB B ( ω , ω ) = (cid:48) (cid:88) P (cid:88) λ λ λ A λ λ B ,λ λ B ,λ λ ω + ω + E λ λ + i [ η + η ] × (cid:20) P λ λ ω + E λ λ + iη − P λ λ ω + E λ λ + iη (cid:21) . (25)We should set η = η in order to turn on all the external fields at the samerate in the adiabatic switching-on procedure. In the imaginary time, τ ≡ it , the time evolution of an operator obeys:ˆ O ( τ ) = e τ H ˆ O e − τ H . We should notice that only in this subsection τ rep-resents an imaginary time. We generalize the periodic property of a bosoniccorrelation function, i.e. f ( τ ) = f ( τ + β ), [27] to a multi-variable correlationfunction as f ( τ , τ , . . . , τ m ). In other words, for an arbitrary choice of τ i , we9ave f ( τ , τ , . . . , τ i , . . . , τ m ) = f ( τ , τ , . . . , τ i + β, . . . , τ m ) . Therefore, we canintroduce the following Fourier and inverse-Fourier transformations: f ( iν , iν , . . . , iν m ) = (cid:90) β dτ (cid:90) β dτ · · · (cid:90) β dτ m f ( τ , τ , . . . , τ m ) e i (cid:80) i ν i τ i ,f ( τ , τ , . . . , τ m ) = (cid:18) β (cid:19) m (cid:88) { ν i } f ( iν , iν , . . . , iν m ) e − i (cid:80) i ν i τ i , (26)where ν i = 0 , ± π/β, ± π/β . . . stands for the bosonic Matsubara frequency(energy). Now, we employ the Matsubara representation for the case of thesecond-order response in order to generalize the first-order Matsubara techniqueto higher-order response functions. The imaginary time-ordered second-orderresponse function is defined by (cid:101) χ AB B ( τ , τ ) = ( − (cid:68) T ˆ B ( − τ ) ˆ B ( − τ ) ˆ A (0) (cid:69) . (27)After performing the time-ordering operation, we arrive at (cid:101) χ AB B ( τ , τ ) = 12 (cid:40) Θ( τ − τ )Θ( − τ ) (cid:68) ˆ B ( − τ ) ˆ B ( − τ ) ˆ A (0) (cid:69) + Θ( τ − τ )Θ( τ ) (cid:68) ˆ A (0) ˆ B ( − τ ) ˆ B ( − τ ) (cid:69) + Θ( τ − τ )Θ( − τ ) (cid:68) ˆ B ( − τ ) ˆ B ( − τ ) ˆ A (0) (cid:69) + Θ( τ − τ )Θ( τ ) (cid:68) ˆ A (0) ˆ B ( − τ ) ˆ B ( − τ ) (cid:69) + Θ( τ )Θ( − τ ) (cid:68) ˆ B ( − τ ) ˆ A (0) ˆ B ( − τ ) (cid:69) + Θ( τ )Θ( − τ ) (cid:68) ˆ B ( − τ ) ˆ A (0) ˆ B ( − τ ) (cid:69) (cid:41) . (28)Since 0 ≤ τ i ≤ β —see Eq. (26)— we have Θ( − τ i ) = 0 and Θ( τ i ) = 1 whichimplies (cid:101) χ AB B ( τ , τ ) = 12 (cid:110) Θ( τ − τ ) (cid:68) ˆ A (0) ˆ B ( − τ ) ˆ B ( − τ ) (cid:69) + Θ( τ − τ ) (cid:68) ˆ A (0) ˆ B ( − τ ) ˆ B ( − τ ) (cid:69) (cid:111) . (29)The above equation can be rewritten in a compact form as follows (cid:101) χ AB B ( τ , τ ) = (cid:48) (cid:88) P Θ( τ − τ ) (cid:68) ˆ A (0) ˆ B ( − τ ) ˆ B ( − τ ) (cid:69) , (30)10here (cid:80) (cid:48)P is the usual intrinsic permutation operation symbol. The thermody-namic average can be taken in the spectral representation: (cid:101) χ AB B ( τ , τ ) = (cid:48) (cid:88) P (cid:88) λ λ λ Θ( τ − τ ) P λ (cid:104) λ | ˆ A (0) | λ (cid:105)(cid:104) λ | ˆ B ( − τ ) | λ (cid:105)(cid:104) λ | ˆ B ( − τ ) | λ (cid:105) . (31)According to the time-evolution convention and considering Eq. (26) for theFourier transformation definition in the imaginary-time domain, we Fouriertransform (cid:101) χ AB B ( τ , τ ) obtaining (cid:101) χ AB B ( iν , iν ) = (cid:48) (cid:88) P (cid:88) λ λ λ P λ A λ λ B ,λ λ B ,λ λ × (cid:90) β dτ (cid:90) β dτ Θ( τ − τ ) e τ E λ λ e τ E λ λ e iν τ e iν τ = (cid:48) (cid:88) P (cid:88) λ λ λ P λ A λ λ B ,λ λ B ,λ λ (cid:40) e βE λ λ [ e βE λ λ − iν + E λ λ )( iν + E λ λ ) − e βE λ λ − iν + E λ λ )( i [ ν + ν ] + E λ λ ) (cid:41) . (32)Using that P λ = e − βE λ / Z , we can simplify the above relationship as follows (cid:101) χ AB B ( iν , iν ) = − (cid:48) (cid:88) P (cid:88) λ λ λ A λ λ B ,λ λ B ,λ λ (cid:40) P λ − P λ ( iν + E λ λ )( iν + E λ λ ) − P λ − P λ ( iν + E λ λ )( i [ ν + ν ] + E λ λ ) (cid:41) . (33)We can write the fraction in the first line of the previous equation in the followingway 1 iν + E λ λ iν + E λ λ = 1 i [ ν + ν ] + E λ λ (cid:20) iν + E λ λ + 1 iν + E λ λ (cid:21) . (34)Using the above identity, we can simplify the expression of (cid:101) χ AB B ( iν , iν ) evenmore and reach (cid:101) χ AB B ( iν , iν ) = − (cid:48) (cid:88) P (cid:88) λ λ λ A λ λ B ,λ λ B ,λ λ i [ ν + ν ] + E λ λ × (cid:26) P λ − P λ iν + E λ λ + P λ − P λ iν + E λ λ − P λ − P λ iν + E λ λ (cid:27) = − (cid:48) (cid:88) P (cid:88) λ λ λ A λ λ B ,λ λ B ,λ λ i [ ν + ν ] + E λ λ (cid:26) P λ − P λ iν + E λ λ − P λ − P λ iν + E λ λ (cid:27) . (35)11fter exchanging λ with λ , we obtain (cid:101) χ AB B ( iν , iν ) = (cid:48) (cid:88) P (cid:88) λ λ λ A λ λ B ,λ λ B ,λ λ i [ ν + ν ] + E λ λ (cid:26) P λ − P λ iν + E λ λ − P λ − P λ iν + E λ λ (cid:27) . (36)By considering the intrinsic permutation operation, we reach the following com-pact form for the spectral representation of the second-order response in theMatsubara frequency domain: (cid:101) χ AB B ( iν , iν ) = (cid:48) (cid:88) P (cid:88) λ λ λ A λ λ B ,λ λ B ,λ λ i [ ν + ν ] + E λ λ (cid:26) P λ λ iν + E λ λ − P λ λ iν + E λ λ (cid:27) . (37)By performing the analytical continuation iν i → ω i + iη i we obtain the physical(retarded) second-order response function given in Eq. (25). Once again, weshould set η = η in order to turn on all the external fields at the same rate inthe adiabatic switching-on procedure. We define the m -th order time-ordered response function as follows: (cid:101) χ AB B ...B m ( τ , τ , . . . , τ m ) = ( − m m ! (cid:68) T ˆ B m ( − τ m ) . . . ˆ B ( − τ ) ˆ B ( − τ ) ˆ A (0) (cid:69) , (38)where τ i denotes an imaginary time. Notice that the ( − m /m ! factor is requiredin order to achieve consistency with the real-time picture, see e.g. Eq. (25). Theintrinsic permutation symmetry is implicitly taken into account through thepresence of the time-ordering operation, T , and the 1 /m ! pre-factor is there toavoid multiple counting. A similar relation for the imaginary-time nonlinearcorrelation function was first reported in Ref. [29]. The response function canbe expressed in terms of a “universal” kernel X ( n ) in the following manner: (cid:101) χ AB ...B n ( iν , . . . , iν n ) = (cid:48) (cid:88) P (cid:88) { λ i } X ( n ) λ ...λ n +1 ( iν , . . . , iν n ) (cid:2) A λ λ n +1 Π ni =1 B i,λ i +1 λ i (cid:3) . (39)The universal kernel, which entirely accounts for the frequency dependence ofthe response, can be evaluated in a recursive manner as described below. In thezero-th order it is given by the statistical occupation factor: X (0) λ = P λ = e − βE λ (cid:80) λ e − βE λ . (40)The first-order case then follows X (1) λ λ ( iν ) = X (0) λ − X (0) λ iν + E λ λ . (41)12ollowing the derivation given above regarding the analytical continuation ofthe second-order response function, we find X (2) λ λ λ ( iν , iν ) = X (1) λ λ ( iν ) − X (1) λ λ ( iν ) i [ ν + ν ] + E λ λ . (42)Similarly, for the third order we obtain X (3) λ λ λ λ ( iν , iν , iν ) = X (2) λ λ λ ( iν , iν ) − X (2) λ λ λ ( iν , iν ) i [ ν + ν + ν ] + E λ λ , (43)and thus eventually for the n -th order we have X ( n ) λ ...λ n +1 ( iν , . . . , iν n ) = X ( n − λ ...λ n ( iν , . . . , iν n − ) − X ( n − λ ...λ n +1 ( iν , . . . , iν n ) i [ ν + · · · + ν n ] + E λ λ n +1 . (44)This is the desired recursion relation. Having the universal X ( n ) response func-tion, one can calculate the n -th order physical response function after properlyincorporating the form-factor part shown in the square bracket in Eq. (39) andsumming on all degrees of freedoms λ i . We have already established that, up tosecond-order, the analytic continuation from time-ordered to causal response iseffected by the replacement, iν i → ω i + iη i . The recursion relation (44) showsthat the same procedure will also work for the n -th order response if it worksfor the n − at all orders . Our proof is considerably simplerthan the one first reported in Ref. [29].
4. Light-matter interaction: gauge transformation
After the above formal analysis of nonlinear response functions, we nowproceed to discuss the issue of gauge invariance in nonlinear response theory,analysing the role of an external electromagnetic field. We then generalize thelinear Ward identity stemming from gauge invariance and charge conservationto higher-order response functions [15, 28].In the interaction representation, the light-matter interaction can be writtenin terms of vector, A ( r , t ), and scalar, Φ( r , t ), potentials in the space-timedomain, ( r , t ):ˆ V ( t ) = − (cid:88) α (cid:90) (cid:18)(cid:90) dλ ˆ J α ( r , t ; λ A ) (cid:19) A α ( r , t ) d r + (cid:90) ˆ n ( r , t )Φ( r , t ) d r , (45)where ˆ J α ( r , t ; A ) is the charge current operator, which, in general, may de-pend on the electromagnetic vector potential A . The integration over the realparameter λ in the first term of Eq. (45) guarantees the fundamental relationbetween the current and the light-matter interaction Hamiltonian:ˆ J α ( r , t ; A ) = δ ˆ V ( t ) δ A α ( r , t ) (46)13or any value of A .Note that the time dependence of the operators ˆ V and ˆ J α originates fromtwo sources: the interaction-picture of the time evolution, see Eq. (2), and theexplicit time dependence of the external potentials. The time dependence ofthe charge-density operator ˆ n stems from the interaction-picture of the timeevolution. The charge-current operator in the interaction picture is given byˆ J α ( r , t ; A ) = ˆ j α ( r , t ) + (cid:88) α (cid:90) d r (cid:48) ˆ κ α α ( r , r (cid:48) , t ) A β ( r (cid:48) , t )+ (cid:88) α α (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) ˆ ξ α α α ( r , r (cid:48) , r (cid:48)(cid:48) , t ) A α ( r (cid:48) , t ) A α ( r (cid:48)(cid:48) , t )+ (cid:88) α α α (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) (cid:90) d r (cid:48)(cid:48)(cid:48) ˆ ζ α α α α ( r , r (cid:48) , r (cid:48)(cid:48) , t ) A α ( r (cid:48) , t ) A α ( r (cid:48)(cid:48) , t ) A α ( r (cid:48)(cid:48)(cid:48) , t )+ . . . , (47)where ˆ H tot = ˆ H + ˆ V . Here, the one-photon current component in the interactionpicture is given byˆ j α ( r , t ) = ˆ J α ( r , t ; A ) | A → = − e it ˆ H δ ˆ H tot δ A α ( r , t ) (cid:12)(cid:12)(cid:12) A → e − it ˆ H . (48)For the two-photon current coupling we haveˆ κ α α ( r , r (cid:48) , t ) = δ ˆ J α ( r , t ; A ) δ A α ( r (cid:48) , t ) (cid:12)(cid:12)(cid:12) A → = − e it ˆ H δ ˆ H tot δ A α ( r , t ) δ A α ( r (cid:48) , t ) (cid:12)(cid:12)(cid:12) A → e − it ˆ H . (49)Similarly, three- and four-photon current couplings read as following:ˆ ξ α α α ( r , r (cid:48) , r (cid:48)(cid:48) , t ) = 12! δ ˆ J α ( r , t ; A ) δ A α ( r (cid:48) , t ) δ A α ( r (cid:48)(cid:48) , t ) (cid:12)(cid:12)(cid:12) A → = − e it ˆ H δ ˆ H tot δ A α ( r , t ) δ A α ( r (cid:48) , t ) δ A α ( r (cid:48)(cid:48) , t ) (cid:12)(cid:12)(cid:12) A → e − it ˆ H (50)andˆ ζ α α α α ( r , r (cid:48) , r (cid:48)(cid:48) , r (cid:48)(cid:48)(cid:48) , t ) = 13! δ ˆ J α ( r , t ; A ) δ A α ( r (cid:48) , t ) δ A α ( r (cid:48)(cid:48) , t ) δ A α ( r (cid:48)(cid:48)(cid:48) , t ) (cid:12)(cid:12)(cid:12) A → = − e it ˆ H δ ˆ H tot δ A α ( r , t ) δ A α ( r (cid:48) , t ) δ A α ( r (cid:48)(cid:48) , t ) δ A α ( r (cid:48)(cid:48)(cid:48) , t ) (cid:12)(cid:12)(cid:12) A → e − it ˆ H . (51)It is common to dub ˆ j α “ paramagnetic ” current operator and ˆ κ α α , ˆ ξ α α α ,and ˆ ζ α α α α “ multi-photon ” current operators (which can have either a dia-magnetic or a paramagnetic contribution to total current depending on the14round-state phase and the order of perturbation; jargon stemming from thetheory of superconductivity [19]). Note, finally, that the time dependence ofˆ j α ( r , t ) (and all multi-photon current operators) stems from the interaction-picture of the time evolution, see Eq. (2).Formally, we have the following perturbative series for the macroscopiccharge density and current N ( r , t ) = (cid:104) ˆ n ( r , t ) (cid:105) + N (1) ( r , t ) + N (2) ( r , t ) + N (3) ( r , t ) + . . . , (52) J (cid:96) ( r , t ) = (cid:68) ˆ j (cid:96) ( r , t ) (cid:69) + J (1) (cid:96) ( r , t ) + J (2) (cid:96) ( r , t ) + J (3) (cid:96) ( r , t ) + . . . . (53)We now apply a gauge transformation to obtain gauge-invariant relations forthe first- and second-order density and current. The gauge transformation isdefined by the following map A (cid:48) ( r , t ) = A ( r , t ) + ∇ Λ( r , t ) , Φ (cid:48) ( r , t ) = Φ( r , t ) − ∂ Λ( r , t ) ∂t , (54)where Λ( r , t ) is an arbitrary smooth function of r and t . Similarly, in Fourierspace: A (cid:48) α ( q , ω ) = A α ( q , ω ) + iq α Λ( q , ω ) , Φ (cid:48) ( q , ω ) = Φ( q , ω ) + iω Λ( q , ω ) . (55)The gauge transformation preserves the electric and magnetic fields, E ( r , t ) = − ∇ Φ( r , t ) − ∂ A ( r , t ) ∂t , B ( r , t ) = ∇ × A ( r , t ) , (56)where we have set c = 1 for the speed of light in vacuum. In Fourier space, thefield components read E α ( q , ω ) = − iq α Φ( q , ω ) + iω A α ( q , ω ) , B α ( q , ω ) = i (cid:88) βγ (cid:15) αβγ q β A γ ( q , ω ) , (57)where (cid:15) αβγ is the Levi-Civita symbol. From now on, we drop the imaginarypart of ω in the notation but keep in mind that ω ≡ ω + i + .According to gauge invariance, the difference between physical observablescalculated in two arbitrary gauges must be zero. This implies: δN ( n ) = N (cid:48) ( n ) − N ( n ) = 0 and δJ ( n ) (cid:96) = J (cid:48) ( n ) (cid:96) − J ( n ) (cid:96) = 0 .
5. Linear response theory: gauge invariance
By truncating the perturbative series up to linear order in the external field(see Appendix B), we reach the following formal relations for the linear charge15ensity and current: N (1) ( q , ω ) = (cid:88) q (cid:48) (cid:40) χ nn ( q , q (cid:48) , ω )Φ( q (cid:48) , ω ) − (cid:88) α χ nj α ( q , q (cid:48) , ω ) A α ( q (cid:48) , ω ) (cid:41) , (58) J (1) (cid:96) ( q , ω ) = (cid:88) q (cid:48) (cid:40) (cid:88) α (cid:104) (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) ) (cid:105) − χ j (cid:96) j α ( q , q (cid:48) , ω ) (cid:105) A α ( q (cid:48) , ω )+ χ j (cid:96) n ( q , q (cid:48) , ω )Φ( q (cid:48) , ω ) (cid:41) . (59)Note that our system is not assumed to be translationally invariant and conse-quently q (cid:48) and q are two independent variables.According to gauge invariance and using Eq. (55), the difference between thephysical charge densities and currents calculated in two different gauges mustbe identically zero: δN (1) ( q , ω ) = i (cid:88) q (cid:48) (cid:40) ωχ nn ( q , q (cid:48) , ω ) − (cid:88) α χ nj α ( q , q (cid:48) , ω ) q (cid:48) α (cid:41) Λ( q (cid:48) , ω ) = 0 ,δJ (1) (cid:96) ( q , ω ) = i (cid:88) q (cid:48) (cid:40) (cid:88) α [ (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) ) (cid:105) − χ j (cid:96) j α ( q , q (cid:48) , ω )] q (cid:48) α + ωχ j (cid:96) n ( q , q (cid:48) , ω ) (cid:41) × Λ( q (cid:48) , ω ) = 0 . (60)In order to fulfill these constraints, the following gauge-invariance identitiesmust be satisfied χ nn ( q , q (cid:48) , ω ) = (cid:88) α q (cid:48) α ω χ nj α ( q , q (cid:48) , ω ) ,χ j (cid:96) n ( q , q (cid:48) , ω ) = (cid:88) α q (cid:48) α ω [ χ j (cid:96) j α ( q , q (cid:48) , ω ) − (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) ) (cid:105) ] . (61)After implementing the above identities in Eq. (58) and using Eq. (57), we canwrite down the following gauge-invariant form of the linear density and current: N (1) ( q , ω ) = − (cid:88) q (cid:48) (cid:88) α χ nj α ( q , q (cid:48) , ω ) E α ( q (cid:48) , ω ) iω , (62) J (1) (cid:96) ( q , ω ) = − (cid:88) q (cid:48) (cid:88) α [ χ j (cid:96) j α ( q , q (cid:48) , ω ) − (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) ) (cid:105) ] E α ( q (cid:48) , ω ) iω . (63)The surprising absence of the magnetic field in these equations is explained bynoting that the finite frequency components of the magnetic field are connectedto the electric field by the two “internal” Maxwell equations, e.g. Faraday’slaw q × E ( q , ω ) = ω B ( q , ω ) and q · B ( q , ω ) = 0, which follow from Eq. (57).16f we expand the conductivity up to linear order in q we can generate termsproportional to the magnetic field.Because of the above relation for the first-order current, we have the followingexpression for the linear conductivity: σ (1) (cid:96)α ( q , q (cid:48) , ω ) = − χ j (cid:96) j α ( q , q (cid:48) , ω ) − (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) ) (cid:105) iω . (64)Utilizing the diamagnetic sum rule for a longitudinal external field E ( q (cid:48) , ω ) || q (cid:48) ,we have (see Ref. [21] and also Appendix C)For longitudinal ˆ j α : (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) ) (cid:105) = lim ω → χ j (cid:96) j α ( q , q (cid:48) , ω ) . (65)For a transverse field E ( q (cid:48) , ω ) ⊥ q (cid:48) the total static current response is finite andtherefore the diamagnetic sum rule is valid only at q (cid:48) = 0:For transverse ˆ j α : (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) = 0) (cid:105) = lim q (cid:48) → lim ω → χ j (cid:96) j α ( q , q (cid:48) , ω ) . (66)The diamagnetic sum rule identity is explicitly proven in Appendix C by usingthe well-known Kubo identity [30, 31] and the continuity relation. This can-cellation clarifies the importance of the diamagnetic component of the currentoperator, i.e. (cid:80) α ˆ κ (cid:96)α A α , in the linear conductivity. Finally, we recall that inthe homogeneous (translational invariant) electron liquid we have q = q (cid:48) . Inthe clean system, the Drude divergence of the dc limit, ω →
0, is achieved inthe reverse order of limits that is first setting q → ω →
0. In thiscase, the paramagnetic response function vanishes, lim ω → lim q → χ j (cid:96) j α →
0, ina perfectly homogeneous electron liquid owing to momentum conservation [21].As a consequence, the diamagnetic contribution leads to the Drude divergenceat zero frequency.The continuity constitution laws for the two- and three-photon diamagneticcurrent couplings read (see Appendix D):[ˆ j β ( r (cid:48) ) , ˆ n ( r )] = (cid:88) α ( − i∂ r α )ˆ κ αβ ( r , r (cid:48) ) , (67)and 12 [[ˆ j γ ( r (cid:48)(cid:48) ) , ˆ n ( r (cid:48) )] , ˆ n ( r )] = (cid:88) αβ ( − i∂ r (cid:48) β )( − i∂ r α ) ˆ ξ αβγ ( r , r (cid:48) , r (cid:48)(cid:48) ) . (68)A similar strategy can be employed to resolve a connection between the equal-time nested commutation of the paramagnetic current and the density operatorwith higher-order n -photon diamagnetic current operators. Let us focus on thetwo-photon diamagnetic current coupling, which, in Fourier space, obeys[ˆ j β ( − q (cid:48) ) , ˆ n ( q )] = (cid:88) α q α ˆ κ αβ ( q , q (cid:48) ) . (69)17t is straightforward to evaluate [ˆ j β ( − q (cid:48) ) , ˆ n ( q )] and obtain (see Appendix D) (cid:88) α q α ˆ κ αβ ( q , q (cid:48) ) = e v β ( − q )ˆ n ( q − q (cid:48) ) . (70)Notice that e > v β ( q ) = ∂ q β H ( q ) isthe velocity operator.
6. Second-order response theory: gauge invariance
In order to obtain the second-order density and current, we need to keep allterms up to quadratic order in the external field. The details of this calculationare available in Appendix B in which we start from the position space andreal time representation. Carrying out lengthy but straightforward calculations(see Appendix B), we reach the following relation for the second-order densityin the Fourier representation: N (2) ( q , ω Σ ) = (cid:88) q q (cid:110) (cid:88) α α Π n ; α α ( Q , Ω ) A α ( q , ω ) A α ( q , ω ) − (cid:88) α Π n ; nα ( Q , Ω )Φ( q , ω ) A α ( q , ω )+ Π n ; nn ( Q , Ω )Φ( q , ω )Φ( q , ω ) (cid:111) , (71)where we have introduced the shortand ( q , q , q , ω , ω ) = ( Q , Ω ). No-tice that Π n ; nn ( Q , Ω ) = (cid:80) (cid:48)P χ nnn ( Q , Ω ), where the Π n ; α α and Π n ; nα response functions are defined byΠ n ; α α ( Q , Ω ) = (cid:48) (cid:88) P (cid:8) χ nj α j α ( Q , Ω ) − χ nκ α α ( Q , Ω ) (cid:9) , Π n ; nα ( Q , Ω ) = χ nnj α ( Q , Ω ) + χ nj α n ( Q , Ω ) . (72)Note that the second-order correlation function χ ABC is defined generically inEq. (13) and (16); see also Eq. (B.8) and Eq. (B.9). Here, (cid:80) (cid:48)P stands for theintrinsic permutation symmetry of each correlation function with respect to thedummy variables [3]. Notice that in our notation “Π A ; B B ... F F . . . ” denotesthe macroscopic value of the “ ˆ A ” operator in response to the “ F i ” externalfields, which are coupled to the “ ˆ B i ” operators.Carrying out the gauge transformation (55) and following similar steps tothose described in the linear-response case, we reach the following relation forthe difference of the second-order densities between the two gauges: δN (2) ( q , ω Σ ) = (cid:88) q (cid:40) i (cid:34)(cid:88) q K ( Q , Ω ; A , Φ) (cid:35) Λ( q , ω ) − (cid:88) q , q K ( Q , Ω )Λ( q , ω )Λ( q , ω ) (cid:41) . (73)18ere, K ( Q , Ω ) = (cid:88) α α Π n ; α α ( Q , Ω ) q ,α q ,α − (cid:88) α Π n ; nα ( Q , Ω ) ω q ,α + Π n ; nn ( Q , Ω ) ω ω (74)and K ( Q , Ω ; A , Φ) = (cid:88) α α (cid:104) Π n ; α α ( Q , Ω ) + Π n ; α α ( Q , Ω ) (cid:105) × q ,α A α ( q , ω ) − (cid:88) α (cid:104) Π n ; nα ( Q , Ω ) ω A α ( q , ω )+ Π n ; nα ( Q , Ω )Φ( q , ω ) q ,α (cid:105) + (cid:104) Π n ; nn ( Q , Ω ) + Π n ; nn ( Q , Ω ) (cid:105) ω Φ( q , ω ) . (75)Notice that in the argument of K we use A and Φ to point out the functionaldependence of K on the external fields. For gauge invariance to hold at second-order perturbation theory, the quantity δN (2) ( q , ω Σ ) must vanish identically.Since δN (2) ( q , ω Σ ) = 0 must be true for any Λ, we need to have (cid:80) q K = 0and K = 0 separately. The condition K = 0 implies (cid:88) α α Π n ; α α ( Q , Ω ) q ,α q ,α − (cid:88) α Π n ; nα ( Q , Ω ) ω q ,α + Π n ; nn ( Q , Ω ) ω ω = 0 . (76)Replacing Eq. (76) in Eq. (71) we obtain a new relation for the second-orderdensity: N (2) ( q , ω Σ ) = (cid:88) q (cid:48) q (cid:48)(cid:48) (cid:88) α α Π n ; α α ( Q , Ω ) E α ( q , ω ) iω E α ( q , ω ) iω + (cid:88) q (cid:48) q (cid:48)(cid:48) (cid:110) (cid:88) α U α ( Q , Ω ) E α ( q , ω ) (cid:111) Φ( q , ω ) iω ω , (77)where U α ( Q , Ω ) = 2 (cid:80) α Π n ; α α ( Q , Ω ) q ,α − Π n ; nα ( Q , Ω ) ω . FromEq. (77), it is evident that there is still a dependence on the scalar potentialthat makes the new representation of N (2) gauge dependent. This is becausewe should also satisfy (cid:80) q K = 0 to reach a fully gauge invariant expressionfor N (2) . However, instead of using (cid:80) q K = 0, we perform another gaugetransformation on Eq. (77): δN (2) ( q , ω Σ ) = − (cid:88) q q (cid:110) (cid:88) α U α ( Q , Ω ) E α ( q , ω ) (cid:111) Λ( q , ω ) ω = 0 . (78)19rom Eq. (78) we conclude that (cid:88) q (cid:88) α U α ( Q , Ω ) E α ( q , ω ) = 0 . (79)By replacing Eq. (79) in Eq. (77), we finally obtain a gauge-invariant expressionfor the second-order density: N (2) ( q , ω Σ ) = (cid:88) q q (cid:88) α α Π n ; α α ( Q , Ω ) E α ( q , ω ) iω E α ( q , ω ) iω . (80)Note that, since Eq. (79) is valid for an arbitrary electric field, it must be U α ( Q , Ω ) = 0 . (81)By plugging Eq. (81) in Eq. (76), we findΠ n ; nn ( Q , Ω ) = (cid:88) α α Π n ; α α ( Q , Ω ) q ,α q ,α ω ω . (82)There is a simpler approach to prove the above identity. Using the scalarpotential gauge, we can write the electric field as E α ( q , ω ) = − iq α Φ( q , ω ).The same electric field is obtained by using a longitudinal vector potential A α ( q , ω ) = − ( q α /ω )Φ( q , ω ). The second-order density fluctuation calculatedin the two gauges must be equal owing to gauge invariance. Therefore, weobtain N (2) ( q , ω Σ ) = Π nnn Φ( q , ω )Φ( q , ω ) = (cid:88) α α Π nα α A α ( q , ω ) A α ( q , ω )= (cid:88) α α Π nα α ( − q ,α Φ( q , ω ) /ω )( − q ,α Φ( q , ω ) /ω ) . (83)Accordingly, we find the same identity given in Eq. (82) after cancelling thescalar potential from the above relation.The above relation is an identity which holds true to ensure the gauge in-variance of the second-order density N (2) ( q , ω Σ ). By using Eqs. (76) and (81)and considering the intrinsic permutation symmetry of Π n ; α α , one can provethat (cid:80) q K = 0 is also fulfilled. This proof is reported in Appendix E.We now proceed to obtain a gauge-invariance relation for the second-ordercurrent. The second-order current in the frequency and wave-vector domain canbe written as follows (see Appendix B) J (2) (cid:96) ( q , ω Σ ) = (cid:88) q q (cid:110) (cid:88) α α Π (cid:96) ; α α ( Q , Ω ) A α ( q , ω ) A α ( q , ω ) − (cid:88) α Π (cid:96) ; nα ( Q , Ω )Φ( q , ω ) A α ( q , ω )+ Π (cid:96) ; nn ( Q , Ω )Φ( q , ω )Φ( q , ω ) (cid:111) , (84)20here Π (cid:96) ; nn ( Q , Ω ) = (cid:80) (cid:48)P χ j α nn ( Q , Ω ) andΠ (cid:96) ; α α ( Q , Ω ) = (cid:48) (cid:88) P (cid:104) χ j (cid:96) j α j α ( Q , Ω ) − χ j (cid:96) κ α α ( q , q , q , ω + ω ) − χ κ (cid:96)α j α ( q , q , q , ω ) + (cid:68) ˆ ξ (cid:96)α α ( q , q , q ) (cid:69) (cid:105) , (85)Π (cid:96) ; nα ( Q , Ω ) = χ j (cid:96) nj α ( Q , Ω ) + χ j (cid:96) j α n ( Q , Ω ) − (cid:48) (cid:88) P χ κ (cid:96)α n ( q , q , q , ω ) . (86)For the following two correlation functions, one needs to be careful in the sym-metrization process, which should be carried out as follows (cid:48) (cid:88) P χ κ (cid:96)α j α ( q , q , q , ω ) = 12 (cid:2) χ κ (cid:96)α j α ( q , q , q , ω ) + χ κ (cid:96)α j α ( q , q , q , ω ) (cid:3) , (87) (cid:48) (cid:88) P χ κ (cid:96)α n ( q , q , q , ω ) = 12 (cid:2) χ κ (cid:96)α n ( q , q , q , ω ) + χ κ (cid:96)α n ( q , q , q , ω ) (cid:3) . (88)After performing the gauge transformation, we reach the following relationfor the second-order gauge-induced current change, δJ (2) (cid:96) ( q , ω Σ ): δJ (2) (cid:96) ( q , ω Σ ) = (cid:88) q i (cid:34)(cid:88) q L ( Q , Ω ; A , Φ) (cid:35) Λ( q , ω ) − (cid:88) q q L ( Q , Ω )Λ( q , ω )Λ( q , ω ) = 0 , (89)where L ( Q , Ω ) = (cid:88) α α Π (cid:96) ; α α ( Q , Ω ) q ,α q ,α − (cid:88) α Π (cid:96) ; nα ( Q , Ω ) ω q ,α + Π (cid:96) ; nn ( Q , Ω ) ω ω (90)and L ( Q , Ω ; A , Φ) = (cid:88) α α (cid:104) Π (cid:96) ; α α ( Q , Ω ) + Π (cid:96) ; α α ( Q , Ω ) (cid:105) A α ( q , ω ) q ,α − (cid:88) α (cid:104) Π (cid:96) ; α n ( Q , Ω ) ω A α ( q , ω )+ Π (cid:96) ; α n ( Q , Ω )Φ( q , ω ) q ,α (cid:105) + (cid:104) Π (cid:96) ; nn ( Q , Ω ) + Π (cid:96) ; nn ( Q , Ω ) (cid:105) Φ( q , ω ) ω . (91)21ecause of gauge invariance, both (cid:80) q L and L must be identically zero, afact that leads to the following relation for L = 0 :Π (cid:96) ; nn ( Q , Ω ) ω ω = (cid:88) α Π (cid:96) ; nα ( Q , Ω ) ω q ,α − (cid:88) α α Π (cid:96) ; α α ( Q , Ω ) q ,α q ,α . (92)Replacing Eq. (92) in Eq. (84) we find J (2) (cid:96) ( q , ω Σ ) = (cid:88) q q (cid:88) α α Π (cid:96) ; α α ( Q , Ω ) E α ( q , ω ) iω E α ( q , ω ) iω + (cid:88) q q (cid:88) α (cid:104) V (cid:96)α ( Q , Ω ) E α ( q , ω ) (cid:105) Φ( q , ω ) iω ω , (93)where V (cid:96)α ( Q , Ω ) = 2 (cid:80) α Π (cid:96) ; α α ( Q , Ω ) q ,α − Π (cid:96) ; nα ( Q , Ω ) ω .Eq. (93) is not a fully gauge-invariant relation for the second-order currentbecause it contains the scalar potential Φ. This is because we have not yet usedthe condition (cid:80) q L = 0. However, for practical reason we instead performanother gauge transformation and later one can show that (cid:80) q L = 0 is alsosatisfied. Accordingly, after performing a gauge transformation to the aboverelation we obtain δJ (2) (cid:96) ( q , ω Σ ) = (cid:88) q q (cid:34)(cid:88) α V (cid:96)α ( Q , Ω ) E α ( q , ω ) (cid:35) Λ( q , ω ) ω = 0 . (94)Because of the gauge invariance property of the physical charge current, weobtain another gauge invariance identity: (cid:88) q (cid:88) α V (cid:96)α ( Q , Ω ) E α ( q , ω ) = 0 . (95)By considering this new identity, we reach the following gauge-invariant relationof the second-order current J (2) (cid:96) ( q , ω Σ ) = (cid:88) q q (cid:88) α α Π (cid:96) ; α α ( Q , Ω ) E α ( q , ω ) iω E α ( q , ω ) iω . (96)Therefore, the second-order conductivity reads as follows: σ (2) (cid:96)α α ( Q , Ω ) ≡ Π (cid:96) ; α α ( Q , Ω ) iω iω . (97)Since Eq. (95) is valid for an arbitrary electric field, we conclude that V (cid:96)α ( Q , Ω ) =0. Using this relation, we can simplify Eq. (92) as followsΠ (cid:96) ; nn ( Q , Ω ) = (cid:88) α α Π (cid:96) ; α α ( Q , Ω ) q ,α q ,α ω ω = i (cid:88) α α q ,α q ,α σ (2) (cid:96)α α ( Q , Ω ) . (98)22 n ˆ n ˆ n ˆ n ˆ n = P `↵ ↵ q ` ! ⌃ q ,↵ q ,↵ ! ! ˆ ⇠ `↵ ↵ + ˆ j ` ˆ j ↵ ˆ j ↵ ˆ j ` ˆ ↵ ↵ ˆ `↵ ˆ j ↵ ˆ n ˆ n ˆ n ˆ n (a) (b) (c) Figure 3: Diagrams for nonlinear response functions in the scalar potential gauge,i.e. Π n ; n...n . (a), (b) and (c) panels correspond to the linear, second-order, and third-orderresponse functions, respectively. Solid lines indicate electronic propagators and solid circlesstand for density vertices. Dashed lines indicate external photon fields. The same result can be obtained in an intuitive way similar to Eq. (83). Theabove equation represents an important Ward identity, which ensures the gaugeinvariance of the second-order current. Similarly to the case of the second-orderdensity, it is easy to prove that (cid:80) q L = 0—this result can be obtained byperforming a calculation similar to that reported in Appendix E for the proofof (cid:80) q K = 0.Finally, we visualize linear, second-order, and third-order response functionsin the scalar and vector potential gauges in Figs. 3 and Fig. 4, respectively. Asit can be seen from Fig. 3, in the scalar gauge, the density response functions ,Π n ; n...n , are given by only one Feynman diagram for any order of perturbation,while in the vector potential gauge, see Fig. 4, we have to evaluate two, four,and eight diagrams for the linear, second-order and third-order current responsefunctions , i.e. Π (cid:96) ; α ...α m , respectively.Formal expression of these diagrams for nonlinear response functions in termsof non-interacting Green’s functions are given in Ref. [14]. Using diagrammaticmethod, nonlinear response functions are explored in two-dimensional linear dis-persive systems known as Dirac materials [13, 32–37] and in strongly-correlatedsystem Ref. [38].
7. Continuity relation and Ward identities
In the absence of external fields there is no current flow (unless the groundstate carries it because of e.g. spontaneous breakdown of time-reversal symme-try). We therefore have ∂ (cid:104) ˆ n ( r , t ) (cid:105) /∂t = 0 and (cid:104) ˆ j α ( r , t ) (cid:105) = 0. In this regard, thecontinuity equation can be written order-by-order in perturbation theory in thefollowing manner: ∇ · J ( m ) ( r , t ) = − ∂N ( m ) ( r , t ) ∂t . (99)In the frequency domain we therefore have: (cid:88) (cid:96) q (cid:96) J ( m ) (cid:96) ( q , ω Σ ) = ω Σ N ( m ) ( q , ω Σ ) . (100)23 n ˆ n ˆ n = P `↵ ↵ q ` ! ⌃ q ,↵ q ,↵ ! ! ˆ ⇠ `↵ ↵ + ˆ j ` ˆ j ↵ ˆ j ↵ ˆ j ` ˆ ↵ ↵ ˆ `↵ ˆ j ↵ ˆ n ˆ n = P `↵ q ` q ↵ ! ˆ `↵ ˆ j ` ˆ j ↵ ˆ j ` ˆ j ↵ ˆ ⇣ `↵ ↵ ↵ ˆ j ` ˆ j ↵ ˆ j ↵ ˆ j ↵ ˆ j ` ˆ j ↵ ˆ ↵ ↵ ˆ j ` ˆ ↵ ↵ ˆ j ↵ ˆ `↵ ˆ j ↵ ˆ j ↵ ˆ j ` ˆ ⇠ ↵ ↵ ↵ ˆ `↵ ˆ ↵ ↵ + + + ˆ ⇠ `↵ ↵ ˆ j ↵ +++ +++ ++ (a)(b)(c) Figure 4: Diagrams for nonlinear response functions in vector potential gauge, i.e. Π (cid:96) ; α ...α m .Solid lines indicate electronic propagators and solid circles stand for the density and currentvertices. (a), (b) and (c) panels correspond to the linear, second-order, and third-order re-sponse function, respectively. Dashed lines indicate external photon fields. The overall signof each bubble obeys the simple rule ( − p − where p is the number of vertices. Notice thatfor all diagrams we must consider the intrinsic permutation symmetrization as discussed inthe main text. Using Eq. (61) in Eq. (100), we have χ nj (cid:96) ( q , q (cid:48) , ω ) = (cid:88) α ( q α /ω ) [ χ j α j (cid:96) ( q , q (cid:48) , ω ) − (cid:104) ˆ κ α(cid:96) ( q , q (cid:48) ) (cid:105) ] . (101)Replacing the above relation in Eq. (61), we arrive at χ nn ( q , q (cid:48) , ω ) = (cid:88) (cid:96)α q (cid:96) q (cid:48) α ω [ χ j (cid:96) j α ( q , q (cid:48) , ω ) − (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) ) (cid:105) ] . (102)We therefore recover the well-known linear-response Ward identity, χ nn ( q , q (cid:48) , ω ) = q · q (cid:48) ω χ J J ( q , q (cid:48) , ω ) , (103)where χ J J ( q , q (cid:48) , ω ) = (cid:88) (cid:96)α q (cid:96) q (cid:48) α q · q (cid:48) [ χ j (cid:96) j α ( q , q (cid:48) , ω ) − (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) ) (cid:105) ] . (104)In a translationally-invariant system q = q (cid:48) and χ J J ( q , q (cid:48) , ω ) represents the longitudinal current-current response function.24imilarly to the linear-response case, we can obtain a “second-order Wardidentity” from the second-order continuity equation. Using Eqs. (80), (96) andEq. (100), we arrive at (cid:88) q q (cid:88) α α (cid:34)(cid:88) (cid:96) q (cid:96) Π (cid:96) ; α α ( Q , Ω ) − ω Σ Π n ; α α ( Q , Ω ) (cid:35) × E α ( q , ω ) iω E α ( q , ω ) iω = 0 . (105)We therefore conclude that Π n ; α α ( Q , Ω ) = (cid:80) (cid:96) q (cid:96) Π (cid:96) ; α α ( Q , Ω ) /ω Σ .Substituting this relation in Eq. (82), we obtain the following second-order Wardidentity:Π n ; nn ( Q , Ω ) = 1 ω Σ (cid:88) (cid:96) q (cid:96) (cid:88) α α q ,α q ,α ω ω Π (cid:96) ; α α ( Q , Ω ) . (106)We can generalize the first- and second-order Ward identities, i.e. Eqs. (102)and (106), to the case of the m -th order response functions as following:Π n ; ...n ( q , q , . . . , q m , ω , . . . , ω m ) = ( − m ω Σ (cid:88) (cid:96) q (cid:96) (cid:88) { α i } q ,α . . . q m,α m ω . . . ω m × Π (cid:96) ; α ...α m ( q , q , . . . , q m , ω , . . . , ω m ) . (107)Notice that in a translationally invariant system, we have q = q Σ = (cid:80) mi =1 q i .Moreover, the m -th order conductivity, the desirable gauge-invariant responsefunction, reads σ ( m ) (cid:96)α ...α m ( q , q , . . . , q m , ω , . . . , ω m ) = Π (cid:96) ; α ...α m ( q , q , . . . , q m , ω , . . . , ω m ) iω . . . iω m . (108)It is useful to discuss the relation between nonlinear conductivities and non-linear density response functions. Considering Eq. (107) and Eq. (108) we findthatΠ n ; ...n ( q , q , . . . , q m , ω , . . . , ω m ) = ( − i ) m ω Σ (cid:88) (cid:96) { α i } q (cid:96) [ q ,α . . . q m,α m ] × σ ( m ) (cid:96)α ...α m ( q , q , . . . , q m , ω , . . . , ω m ) . (109)This relation provides a gauge-covariant prescription in order to relate the dy-namical non-local conductivity to the density response function at each orderof perturbation theory. In the optical (electric-dipole) approximation we canneglect the wave-vector dependence of the conductivity and therefore the wave-25ector expansion of the density response function in the electric-dipole approx-imation reads as following:Π n ; ...n ( q , . . . , q m , ω , . . . , ω m ) ≈ ( − i ) m ω Σ (cid:88) (cid:96) { α i } q Σ ,(cid:96) [ q ,α . . . q m,α m ] × σ ( m ) (cid:96)α ...α m ( ω , . . . , ω m ) . (110)In a translationally-invariant system we have q = q Σ . The latter is therefore not an independent variable, as indicated by the notation “ q Σ ,(cid:96) [ q ,α . . . q m,α m ]”.The quantity σ ( m ) (cid:96)α ...α m ( ω , . . . , ω m ), which is a shorthand for σ ( m ) (cid:96)α ...α m | { q i → } ,is the m -th order optical conductivity.We can go beyond the electric-dipole approximation and consider electric-quadrupole and magnetic dipole contributions by retaining the wave-vector de-pendence of the conductivity up to linear order:Π n ; ...n ( q , . . . , q m , ω , . . . , ω m ) ≈ ( − i ) m ω Σ (cid:88) (cid:96) { α i } q (cid:96) [ q ,α . . . q m,α m ] × (cid:40) σ ( m ) (cid:96)α ...α m ( ω , . . . , ω m ) + (cid:88) i (cid:88) α i q i,α i ∂σ ( m ) (cid:96)α ...α m ( q , . . . , q m , ω , . . . , ω m ) ∂q i,α i (cid:12)(cid:12)(cid:12) { q i → } (cid:41) . (111)The above relation was used [13] in studying nonlinear plasmonic effects ingraphene, which is a centro-symmetric material. Because of the latter property,higher multipole contributions are important in this material [13].
8. Summary
In this Article, we have used the equilibrium ˆ S -matrix approach in orderto set up a theory that allows the calculation of nonlinear response functionsin both real- and imaginary-time domains. We have discussed an analyticalcontinuation procedure for the nonlinear response functions, which provides aprescription in order to obtain the retarded response function from the Matsub-ara one.A large fraction of this work was devoted to the analysis of gauge invari-ance in the context of the nonlinear response theory. An inhomogeneous vectorpotential gauge is a complete gauge to take into account contributions from allmultipoles, such as electric-dipole, electric-quadruple, and magnetic-dipole con-tributions. As a result of gauge invariance, a set of nonlinear Ward identitiesare obtained. Acknowledgment
This work was partially supported by the European Union’s Horizon 2020research and innovation programme under grant No. 881603 (GrapheneCore3).26he work of M. I. K. is supported by European Research Council via SynergyGrant No. 854843 - FASTCORR. G. V. was supported by the U.S. Depart-ment of Energy (Office of Science) under Grant DE- FG02-05ER46203. H. R.acknowledges the support from the Swedish Research Council (VR 2018-04252).
References [1] L. D. Landau and E. M. Lifshitz,
Course of Theoretical Physics: Electro-dynamics of Continuos Media (Pergamon, New York, 1984).[2] Y. R. Shen,
The Principles of Nonlinear Optics (Wiley, New York, 1984)[3] P. N. Butcher and D. Cotter,
The Elements of Nonlinear Optics (CambridgeUniversity Press, Cambridge, 1990).[4] R. W. Boyd,
Nonlinear Optics (Academic Press, Cambridge, 2008).[5] S. Naguleswaran and G. E. Stedman, J. Phys. B: At. Mol. Opt. Phys. ,4027 (1996).[6] A. D. Bandrauk, F. Fillion-Gourdeau, and E. Lorin, J. Phys. B: At. Mol.Opt. Phys. , 153001 (2013).[7] V. N. Genkin and P. M. Mednis, Sov. Phys. JETP , 609 (1968).[8] C. Aversa and J. E. Sipe, Phys. Rev. B , 14636 (1995).[9] G. B. Ventura, D. J. Passos, J. M. B. Lopes dos Santos, J. M. Viana ParenteLopes, and N. M. R. Peres, Phys. Rev. B , 035431 (2017).[10] A. S¨ayn¨atjoki, L. Karvonen, H. Rostami, A. Autere, S. Mehravar, A. Lom-bardo, R. A. Norwood, T. Hasan, N. Peyghambarian, H. Lipsanen, K. Kieu,A. C. Ferrari, M. Polini, and Z. Sun, Nat. Commun. , 893 (2017).[11] S. A. Mikhailov, Phys. Rev. B , 085403 (2016).[12] H. Rostami and M. Polini, Phys. Rev. B , 161411(R) (2016).[13] H. Rostami, M. I. Katsnelson, and M. Polini, Phys. Rev. B , 035416(2017).[14] D. E. Parker, T. Morimoto, J. Orenstein, and J. E. Moore, Phys. Rev. B , 045121 (2019).[15] J.L. Cheng, N. Vermeulen, and J. E. Sipe, Sci Rep. , 43843 (2017).[16] P. Apell, Phys. Scr. , 211 (1983).[17] V. Y. Chernyak, P. Saurabh, and S. Mukamel, J. Chem. Phys. , 164107(2015). 2718] W. E. Lamb, Jr., R. R. Schlicher, and M. O. Scully, Phys. Rev. A , 2763(1987).[19] J. R. Schrieffer, Theory of Superconductivity , (W.A. Benjamin, New York,1964)[20] D. Pines and P. Nozi´eres,
The Theory of Quantum Liquids (W.A. Benjamin,Inc., New York, 1966).[21] G. F. Giuliani and G. Vignale,
Quantum Theory of the Electron Liquid (Cambridge University Press, Cambridge, 2005).[22] A. A. Grinberg and S. Luryi, Phys. Rev. B , 87 (1988).[23] M. M. Glazov and S. D. Ganichev, Phys. Rep. ,101-138 (2014).[24] Lundeberg, Mark B. and Gao, Yuanda and Asgari, Reza and Tan, Chengand Van Duppen, Ben and Autore, Marta and Alonso-Gonz´alez, Pabloand Woessner, Achim and Watanabe, Kenji and Taniguchi, Takashi andHillenbrand, Rainer and Hone, James and Polini, Marco and Koppens,Frank H. L. Science , 187 (2017).[25] F. A. Hopf, P. Meystre, M. O. Scully, and W. H. Louisell, Opt. Commun. , 413 (1976).[26] B. McNeil and N. Thompson, Nat. Photonics , 814 (2010).[27] G. D. Mahan, Many-Particle Physics (Springer US, New York, 1981).[28] M. Peskin and D. V. Schroeder,
An Introduction to Quantum Field Theory ,(Avalon Publishing,New York, 1995).[29] W. A. B. Evans, Proceedings of the Physical Society , 723 (1966).[30] R. C. Clark and G. H. Derrick, Mathematical Methods in Solid State andSuperfluid Theory , (Springer US, New York, 1968).[31] S. G. Louie and M. L. Cohen,
Conceptual Foundations of Materials: AStandard Model for Ground- and Excited-State Properties (Elsevier, NewYork, 2006).[32] T. O. Wehling, A. Huber, A. I. Lichtenstein, and M. I. Katsnelson, Phys.Rev. B , 041404(R) (2015).[33] M. Vandelli, M. I. Katsnelson, and E. A. Stepanov, Phys. Rev. B , 165432(2019).[34] E. A. Stepanov, S. V. Semin, C. R. Woods, M. Vandelli, A. V. Kimel, K.S. Novoselov, and M. I. Katsnelson ACS Appl. Mater. Interfaces , 27758(2020).[35] H. Rostami and V. Juriˇci´c, Phys. Rev. Research , 013069 (2020).2836] H. Rostami and E. Cappelluti, arXiv:2007.08282 (2020).[37] H. Rostami and E. Cappelluti, arXiv:2011.03824 (2020).[38] M. I. Katsnelson and A. I. Lichtenstein, J. Phys.: Condens. Matter ,382201 (2010). Appendix A. Perturbation theory
The macroscopic value of a generic operator ˆ O ( t ) reads as following: O ( t ) = (cid:104) ψ ( t ) | ˆ O ( t ) | ψ ( t ) (cid:105) = (cid:104) ψ | ˆ S † ( t, −∞ ) ˆ O ( t ) ˆ S ( t, −∞ ) | ψ (cid:105) . (A.1)Thanks to the Dyson expansion, we haveˆ S ( t, −∞ ) = 1 + ∞ (cid:88) m =1 ( − i ) m D m ˆ V ( t ) ˆ V ( t ) . . . ˆ V ( t m )ˆ S † ( t, −∞ ) = 1 + ∞ (cid:88) m =1 i m D m ˆ V ( t m ) ˆ V ( t m − ) . . . ˆ V ( t ) , (A.2)where we define D m as follows D m ≡ (cid:90) t −∞ dt (cid:90) t −∞ dt · · · (cid:90) t m − −∞ dt m . (A.3)Using Eq. (A.1), the macroscopic value of ˆ O ( t ) is given by O ( t ) = (cid:68) (cid:110) i D ˆ V ( t ) + i D ˆ V ( t ) ˆ V ( t ) + i D ˆ V ( t ) ˆ V ( t ) ˆ V ( t ) + . . . (cid:111) ˆ O ( t ) (cid:110) − i D ˆ V ( t ) + i D ˆ V ( t ) ˆ V ( t ) − i D ˆ V ( t ) ˆ V ( t ) ˆ V ( t ) + . . . (cid:111) (cid:69) = (cid:68) ˆ O ( t ) (cid:69) + i D (cid:68)(cid:104) ˆ V ( t ) , ˆ O ( t ) (cid:105)(cid:69) + i (cid:68) D ˆ V ( t ) ˆ V ( t ) ˆ O ( t )+ ˆ O ( t ) D ˆ V ( t ) ˆ V ( t ) − D ˆ V ( t ) ˆ O ( t ) D ˆ V ( t ) (cid:69) + i (cid:68) D ˆ V ( t ) ˆ V ( t ) ˆ V ( t ) ˆ O ( t ) − ˆ O ( t ) D ˆ V ( t ) ˆ V ( t ) ˆ V ( t )+ D ˆ V ( t ) ˆ O ( t ) D ˆ V ( t ) ˆ V ( t ) − D ˆ V ( t ) ˆ V ( t ) ˆ O ( t ) D ˆ V ( t ) (cid:69) + . . . . (A.4)In order to simplify the terms proportional to i in the above relation, we needto rewrite D ˆ V ( t ) ˆ O ( t ) D ˆ V ( t ) in a proper form. To do so, we divide the squaredomain of the integral into two triangular domains (see Fig. A.5). We therefore29 ∞ ⋯ − ∞ ⋯ $$ $ % $ & = +−∞ ⋯ − ∞ ⋯ $$ $ % $ & −∞ ⋯ − ∞ ⋯ $$ $ % $ & Figure A.5: From the square domain into two triangle ones. find D ˆ V ( t ) ˆ O ( t ) D ˆ V ( t ) = (cid:90) t −∞ dt (cid:90) t −∞ dt ˆ V ( t ) ˆ O ( t ) ˆ V ( t )= (cid:90) t −∞ dt (cid:90) t −∞ dt ˆ V ( t ) ˆ O ( t ) ˆ V ( t )+ (cid:90) t −∞ dt (cid:90) t −∞ dt ˆ V ( t ) ˆ O ( t ) ˆ V ( t )= (cid:90) t −∞ dt (cid:90) t −∞ dt (cid:110) ˆ V ( t ) ˆ O ( t ) ˆ V ( t ) + ˆ V ( t ) ˆ O ( t ) ˆ V ( t ) (cid:111) . (A.5)Using the above identity, we reach the following compact form for the second-order perturbative terms in Eq. (A.4): D ˆ V ( t ) ˆ V ( t ) ˆ O ( t ) + ˆ O ( t ) D ˆ V ( t ) ˆ V ( t ) − D ˆ V ( t ) ˆ O ( t ) D ˆ V ( t )= (cid:90) t −∞ dt (cid:90) t −∞ dt (cid:104) ˆ V ( t ) , (cid:104) ˆ V ( t ) , ˆ O ( t ) (cid:105)(cid:105) . (A.6)Similarly, one can find the following simplified expression for the third-orderperturbation D ˆ V ( t ) ˆ V ( t ) ˆ V ( t ) O ( t ) − O ( t ) D ˆ V ( t ) ˆ V ( t ) ˆ V ( t ) + D ˆ V ( t ) O ( t ) D ˆ V ( t ) ˆ V ( t ) − D ˆ V ( t ) ˆ V ( t ) O ( t ) D ˆ V ( t )= (cid:90) t −∞ dt (cid:90) t −∞ dt (cid:90) t −∞ dt (cid:104) ˆ V ( t ) , (cid:104) ˆ V ( t ) , (cid:104) ˆ V ( t ) , O ( t ) (cid:105)(cid:105)(cid:105) . (A.7)30ventually, after plugging Eqs. (A.6) and (A.7) in Eq. (A.4), we arrive at thefollowing relation for O ( t ): O ( t ) = (cid:68) ˆ O ( t ) (cid:69) + i (cid:90) t −∞ dt (cid:68)(cid:104) ˆ V ( t ) , ˆ O ( t ) (cid:105)(cid:69) + i (cid:90) t −∞ dt (cid:90) t −∞ dt (cid:68)(cid:104) ˆ V ( t ) , (cid:104) ˆ V ( t ) , ˆ O ( t ) (cid:105)(cid:105)(cid:69) + i (cid:90) t −∞ dt (cid:90) t −∞ dt (cid:90) t −∞ dt (cid:68)(cid:104) ˆ V ( t ) , (cid:104) ˆ V ( t ) , (cid:104) ˆ V ( t ) , ˆ O ( t ) (cid:105)(cid:105)(cid:105)(cid:69) + . . . . (A.8)Introducing the new variables τ i = t − t i , we have O ( t ) = (cid:68) ˆ O ( t ) (cid:69) + i (cid:90) ∞ dτ (cid:68)(cid:104) ˆ V ( t − τ ) , ˆ O ( t ) (cid:105)(cid:69) + i (cid:90) ∞ dτ (cid:90) ∞ τ dτ (cid:68)(cid:104) ˆ V ( t − τ ) , (cid:104) ˆ V ( t − τ ) , ˆ O ( t ) (cid:105)(cid:105)(cid:69) + i (cid:90) ∞ dτ (cid:90) ∞ τ dτ (cid:90) ∞ τ dτ (cid:68)(cid:104) ˆ V ( t − τ ) , (cid:104) ˆ V ( t − τ ) , (cid:104) ˆ V ( t − τ ) , ˆ O ( t ) (cid:105)(cid:105)(cid:105)(cid:69) + . . . . (A.9)Using the Heaviside function Θ( τ ), we reach O ( t ) = (cid:68) ˆ O ( t ) (cid:69) + i (cid:90) ∞−∞ dτ Θ( τ ) (cid:68)(cid:104) ˆ V ( t − τ ) , ˆ O ( t ) (cid:105)(cid:69) + i (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ Θ( τ − τ )Θ( τ ) (cid:68)(cid:104) ˆ V ( t − τ ) , (cid:104) ˆ V ( t − τ ) , ˆ O ( t ) (cid:105)(cid:105)(cid:69) + i (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ Θ( τ − τ )Θ( τ − τ )Θ( τ ) × (cid:68)(cid:104) ˆ V ( t − τ ) , (cid:104) ˆ V ( t − τ ) , (cid:104) ˆ V ( t − τ ) , ˆ O ( t ) (cid:105)(cid:105)(cid:105)(cid:69) + . . . . (A.10) Appendix B. Position and real-time representation of the linear andsecond-order responses
First of all, we calculate the following terms which contain different pertur-bative contributions. The average of the current operator reads (cid:68) ˆ J (cid:96) ( r , t ; A ) (cid:69) = (cid:68) ˆ j (cid:96) ( r , t ) (cid:69) + (cid:90) d r (cid:48) (cid:88) α (cid:104) ˆ κ (cid:96)α ( r , r (cid:48) , t ) (cid:105) A α ( r (cid:48) , t )+ (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) (cid:88) α α (cid:68) ˆ ξ (cid:96)α α ( r , r , r (cid:48)(cid:48) , t ) (cid:69) A α ( r (cid:48) , t ) A α ( r (cid:48)(cid:48) , t ) + . . . . (B.1)31n the canonical ensemble we have [ˆ ρ, ˆ H ] = 0 where ˆ ρ is the density matrix.Therefore, we can show that (cid:104) ˆ κ (cid:96)α ( r , r (cid:48) , t ) (cid:105) = Tr (cid:104) ˆ ρe it ˆ H ˆ κ (cid:96)α ( r , r (cid:48) ) e − it ˆ H (cid:105) = Tr [ˆ ρ ˆ κ (cid:96)α ( r , r (cid:48) )] = (cid:104) ˆ κ (cid:96)α ( r , r (cid:48) ) (cid:105) . (B.2)and similarly, (cid:104) ˆ ξ (cid:96)α α ( r , r (cid:48) , r (cid:48)(cid:48) , t ) (cid:105) = (cid:104) ˆ ξ (cid:96)α α ( r , r (cid:48) , r (cid:48)(cid:48) ) (cid:105) . Eventually, we have (cid:68) ˆ J (cid:96) ( r , t ; A ) (cid:69) = (cid:68) ˆ j (cid:96) ( r , t ) (cid:69) + (cid:90) d r (cid:48) (cid:88) α (cid:104) ˆ κ (cid:96)α ( r , r (cid:48) ) (cid:105) A α ( r (cid:48) , t )+ (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) (cid:88) α α (cid:68) ˆ ξ (cid:96)α α ( r , r , r (cid:48)(cid:48) ) (cid:69) A α ( r (cid:48) , t ) A α ( r (cid:48)(cid:48) , t ) + . . . . (B.3)The other required terms are the following commutators (cid:68)(cid:104) ˆ V ( t ) , ˆ J (cid:96) ( r , t ; A ) (cid:105)(cid:69) = (cid:90) d r (cid:40) (cid:68)(cid:104) ˆ n ( r , t ) , ˆ J (cid:96) ( r , t ; A ) (cid:105)(cid:69) Φ( r , t ) − (cid:88) α (cid:68)(cid:104) ˆ J α ( r , t ; A ) , ˆ J (cid:96) ( r , t ; A ) (cid:105)(cid:69) A α ( r , t ) (cid:41) = (cid:90) d r (cid:40) (cid:68)(cid:104) ˆ n ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) Φ( r , t ) − (cid:88) α (cid:68)(cid:104) ˆ j α ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) A α ( r , t )+ (cid:90) d r (cid:48) (cid:88) α (cid:104) [ˆ n ( r , t ) , ˆ κ (cid:96)α ( r , r (cid:48) , t )] (cid:105) Φ( r , t ) A α ( r (cid:48) , t ) − (cid:90) d r (cid:48) (cid:88) α α (cid:68)(cid:104) ˆ j α ( r , t ) , ˆ κ (cid:96)α ( r , r (cid:48) , t ) (cid:105)(cid:69) A α ( r , t ) A α ( r (cid:48) , t ) − (cid:90) d r (cid:48) (cid:88) α α (cid:68)(cid:104) ˆ κ α α ( r , r (cid:48) , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) × A α ( r , t ) A α ( r (cid:48) , t ) (cid:41) + . . . , (B.4)32n which ˆ V ( t ) is defined in Eq. (45) as the light-matter interaction. Moreover,we calculate the following necessary term (cid:68)(cid:104) ˆ V ( t ) , (cid:104) ˆ V ( t ) , ˆ J (cid:96) ( r , t ; A ) (cid:105)(cid:105)(cid:69) = (cid:90) d r (cid:90) d r (cid:40)(cid:68)(cid:104) ˆ n ( r , t ) , (cid:104) ˆ n ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:105)(cid:69) Φ( r , t )Φ( r , t ) − (cid:88) α (cid:68)(cid:104) ˆ j α ( r , t ) , (cid:104) ˆ n ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:105)(cid:69) Φ( r , t ) A α ( r , t ) − (cid:88) α (cid:68)(cid:104) ˆ n ( r , t ) , (cid:104) ˆ j α ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:105)(cid:69) A α ( r , t )Φ( r , t )+ (cid:88) α α (cid:68)(cid:104) ˆ j α ( r , t ) , (cid:104) ˆ j α ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:105)(cid:69) × A α ( r , t ) A α ( r , t ) (cid:41) + . . . . (B.5)In a similar way we can write down the corresponding expressions for (cid:68)(cid:104) ˆ V ( t ) , ˆ n ( r , t ) (cid:105)(cid:69) and (cid:68)(cid:104) ˆ V ( t ) , (cid:104) ˆ V ( t ) , ˆ n ( r , t ) (cid:105)(cid:105)(cid:69) , which are necessary to obtain the first- andsecond-order densities. Appendix B.1. First- and second-order currents
By plugging Eqs. (B.3), (B.4), and (B.5) in Eq. (A.10), we obtain the fol-lowing expression for the first-order current: J (1) (cid:96) ( r , t ) = (cid:90) d r (cid:48) (cid:88) α (cid:104) ˆ κ (cid:96)α ( r , r (cid:48) ) (cid:105) A α ( r (cid:48) , t )+ i (cid:90) ∞−∞ dτ (cid:90) d r Θ( τ ) (cid:40) (cid:68)(cid:104) ˆ n ( r , t − τ ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) Φ( r , t − τ ) − (cid:88) α (cid:68)(cid:104) ˆ j α ( r , t − τ ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) A α ( r , t − τ ) (cid:41) . (B.6)Using the linear-response definition given in Eq. (12), we rewrite the previousrelation as following: J (1) (cid:96) ( r , t ) = (cid:90) d r (cid:48) (cid:40) (cid:88) α (cid:104) (cid:104) ˆ κ (cid:96)α ( r , r (cid:48) ) (cid:105) A α ( r (cid:48) , t ) − (cid:90) ∞−∞ dτ χ j (cid:96) j α ( r , r (cid:48) , τ ) A α ( r (cid:48) , t − τ ) (cid:105) + (cid:90) ∞−∞ dτ χ j (cid:96) n ( r , r (cid:48) , τ )Φ( r , t − τ ) (cid:41) . (B.7)33n a similar way, the second-order current in position space and time domaincan be written as follows J (2) (cid:96) ( r , t ) = (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) (cid:88) α α (cid:68) ˆ ξ (cid:96)α α ( r , r (cid:48) , r (cid:48)(cid:48) ) (cid:69) A α ( r (cid:48) , t ) A α ( r (cid:48)(cid:48) , t )+ i (cid:90) ∞−∞ dτ (cid:90) d r (cid:90) d r (cid:48) Θ( τ ) (cid:40)(cid:88) α (cid:104) [ˆ n ( r , t − τ ) , ˆ κ (cid:96)α ( r , r (cid:48) , t )] (cid:105) Φ( r , t − τ ) A α ( r (cid:48) , t ) − (cid:88) α α (cid:68)(cid:104) ˆ j α ( r , t − τ ) , ˆ κ (cid:96)α ( r , r (cid:48) , t ) (cid:105)(cid:69) A α ( r , t − τ ) A α ( r (cid:48) , t ) − (cid:88) α α (cid:68)(cid:104) ˆ κ α α ( r , r (cid:48) , t − τ ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) A α ( r , t − τ ) A α ( r (cid:48) , t − τ ) (cid:41) + i (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ (cid:90) d r (cid:90) d r Θ( τ − τ )Θ( τ ) (cid:40)(cid:104) ˆ n ( r , t − τ ) , (cid:104) ˆ n ( r , t − τ ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:105) Φ( r , t − τ )Φ( r , t − τ ) − (cid:88) α (cid:104) ˆ j α ( r , t − τ ) , (cid:104) ˆ n ( r , t − τ ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:105) Φ( r , t − τ ) A α ( r , t − τ ) − (cid:88) α (cid:104) ˆ n ( r , t − τ ) , (cid:104) ˆ j α ( r , t − τ ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:105) A α ( r , t − τ )Φ( r , t − τ )+ (cid:88) α α (cid:104) ˆ j α ( r , t − τ ) , (cid:104) ˆ j α ( r , t − τ ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:105) A α ( r , t − τ ) A α ( r , t − τ ) (cid:41) . (B.8)34sing Eqs. (12) and (13), we rewrite the previous relation as follows J (2) (cid:96) ( r , t ) = (cid:90) d r (cid:90) d r (cid:40) (cid:88) α α (cid:68) ˆ ξ (cid:96)α α ( r , r , r ) (cid:69) A α ( r , t ) A α ( r , t )+ (cid:90) ∞−∞ dτ (cid:34) (cid:88) α χ κ (cid:96)α n ( r , r , r , τ ) A α ( r , t )Φ( r , t − τ ) − (cid:88) α α χ κ (cid:96)α j α ( r , r , r , τ ) A α ( r , t − τ ) A α ( r , t ) − (cid:88) α α χ j (cid:96) κ α α ( r , r , r , τ ) A α ( r , t − τ ) A α ( r , t − τ ) (cid:35) + (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ (cid:34) χ j (cid:96) nn ( r , r , r , τ , τ )Φ( r , t − τ )Φ( r , t − τ ) − (cid:88) α χ j (cid:96) nj α ( r , r , r , τ , τ )Φ( r , t − τ ) A α ( r , t − τ ) − (cid:88) α χ j (cid:96) j α n ( r , r , r , τ , τ ) A α ( r , t − τ )Φ( r , t − τ )+ (cid:88) α α χ j (cid:96) j α j α ( r , r , r , τ , τ ) A α ( r , t − τ ) A α ( r , t − τ ) (cid:35)(cid:41) . (B.9)By performing the Fourier transformation with respect to time time— Eqs. (15)and (16)—and position as explained in Appendix B.3, we obtain the finalexpressions for the first- and second-order currents presented in the main text. Appendix B.2. First- and second-order densities
The first-order density is given by: N (1) ( r , t ) = i (cid:90) ∞−∞ dτ (cid:90) d r (cid:48) Θ( τ ) (cid:104) [ˆ n ( r (cid:48) , t − τ ) , ˆ n ( r , t )] (cid:105) Φ( r (cid:48) , t − τ ) − i (cid:88) α (cid:90) ∞−∞ dτ (cid:90) d r (cid:48) Θ( τ ) (cid:68)(cid:104) ˆ j α ( r (cid:48) , t − τ ) , ˆ n ( r , t ) (cid:105)(cid:69) A α ( r (cid:48) , t − τ ) . (B.10)Using the linear-response definition given in Eq. (12), we rewrite the previousrelation as follows N (1) ( r , t ) = (cid:90) ∞−∞ dτ (cid:90) d r (cid:48) χ nn ( r , r (cid:48) , τ )Φ( r (cid:48) , t − τ ) − (cid:88) α (cid:90) ∞−∞ dτ (cid:90) d r (cid:48) χ nj α ( r , r (cid:48) , τ ) A α ( r (cid:48) , t − τ ) . (B.11)35he second-order density is given by N (2) ( r , t ) = − i (cid:88) α α (cid:90) ∞−∞ dτ (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) Θ( τ ) (cid:104) [ˆ κ α α ( r (cid:48) , r (cid:48)(cid:48) , t − τ ) , ˆ n ( r , t )] (cid:105)× A α ( r (cid:48) , t − τ ) A α ( r (cid:48)(cid:48) , t − τ )+ i (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) Θ( τ − τ )Θ( τ ) (cid:40) (cid:104) [ˆ n ( r (cid:48)(cid:48) , t − τ ) , [ˆ n ( r (cid:48) , t − τ ) , ˆ n ( r , t )]] (cid:105) Φ( r (cid:48) , t − τ )Φ( r (cid:48)(cid:48) , t − τ ) − (cid:88) α (cid:68)(cid:104) ˆ n ( r (cid:48)(cid:48) , t − τ ) , (cid:104) ˆ j α ( r (cid:48) , t − τ ) , ˆ n ( r , t ) (cid:105)(cid:105)(cid:69) A α ( r (cid:48) , t − τ )Φ( r (cid:48)(cid:48) , t − τ ) − (cid:88) α (cid:68)(cid:104) ˆ j α ( r (cid:48)(cid:48) , t − τ ) , [ˆ n ( r (cid:48) , t − τ ) , ˆ n ( r , t )] (cid:105)(cid:69) Φ( r (cid:48) , t − τ ) A α ( r (cid:48)(cid:48) , t − τ )+ (cid:88) α α (cid:68)(cid:104) ˆ j α ( r (cid:48)(cid:48) , t − τ ) , (cid:104) ˆ j α ( r (cid:48) , t − τ ) , ˆ n ( r , t ) (cid:105)(cid:105)(cid:69) × A α ( r (cid:48) , t − τ ) A α ( r (cid:48)(cid:48) , t − τ ) (cid:41) . (B.12)By using the definitions given in Eqs. (12) and (13), we arrive at N (2) ( r , t ) = (cid:90) d r (cid:90) d r (cid:34) − (cid:90) ∞−∞ dτ (cid:88) α α χ nκ α α ( r , r , τ ) A α ( r , t − τ ) A α ( r , t − τ )+ (cid:90) ∞−∞ dτ (cid:90) ∞−∞ dτ (cid:40) χ nnn ( r , r , r , τ , τ )Φ( r , t − τ )Φ( r , t − τ ) − (cid:88) α χ nj α n ( r , r , r , τ , τ ) A α ( r , t − τ )Φ( r , t − τ ) − (cid:88) α χ nnj α ( r , r , r , τ , τ )Φ( r , t − τ ) A α ( r , t − τ )+ (cid:88) α α χ nj α j α ( r , r , r , τ , τ ) A α ( r , t − τ ) A α ( r , t − τ ) (cid:41)(cid:35) . (B.13)By performing the Fourier transformation with respect to time, Eqs. (15) and(16), and position as explained in Appendix B.3, we obtain the final expressionsfor the first- and second-order densities presented in the main text. Appendix B.3. Fourier transformation
In order to go from the position space r to the wave-vector space q , weperform a Fourier transformation. For a general expression of the type G ( r ) = (cid:88) { r j } K ( r , r , . . . r m ) F ( r ) . . . F m ( r m ) , (B.14)36he Fourier-transformed counterpart reads as following: G ( q ) = (cid:88) { q j } K ( q , q , . . . q m ) F ( q ) . . . F m ( q m ) , (B.15)where G ( q ) = (cid:88) r G ( r ) e − i q · r , F j ( q j ) = (cid:88) r j F j ( q j ) e − i q j · r j , K ( q , q , . . . q m ) = (cid:88) r , { r j } K ( r , r , . . . r m ) exp (cid:16) − i q · r + i m (cid:88) j q j · r j (cid:17) . (B.16)We notice that for the case of translationally-invariant system we have q = q Σ = (cid:80) mj q j . Appendix C. Diamagnetic contribution to the linear conductivity:Proof for Eq. (66)
In the scalar potential gauge, we have the following relation for the first-ordercurrent: J (1) (cid:96) ( r , t ) = i (cid:90) t −∞ dt (cid:90) d r (cid:68)(cid:104) ˆ n ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) Φ( r , t ) . (C.1)By using the cyclic property of the trace, we can prove that: (cid:68)(cid:104) ˆ n ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) = Tr (cid:104) ˆ ρ (cid:104) ˆ n ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:105) = Tr (cid:104) [ˆ ρ, ˆ n ( r , t )] ˆ j (cid:96) ( r , t ) (cid:105) = (cid:90) β dλ (cid:68)(cid:104) ˆ n ( r , t − iλ ) , ˆ H (cid:105) ˆ j (cid:96) ( r , t ) (cid:69) . (C.2)Notice that in the last equality of the above relation, we have used the Kuboidentity [30, 31]: (cid:104) ˆ ρ, ˆ O ( t ) (cid:105) = ˆ ρ (cid:90) β dλ (cid:104) ˆ O ( t − iλ ) , ˆ H (cid:105) , (C.3)where ˆ O ( t − iλ ) = e λ ˆ H ˆ O ( t ) e − λ ˆ H . (C.4)By performing a straightforward calculation one can show that (cid:68)(cid:104) ˆ n ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) = (cid:90) β dλ (cid:68)(cid:104) ˆ n ( r , , ˆ H (cid:105) ˆ j (cid:96) ( r , t − t + iλ ) (cid:69) . (C.5)According to the continuity equation, we have d ˆ n ( r , dt = − i [ˆ n ( r , , ˆ H ] = − ∇ · ˆ j ( r , . (C.6)37t is important to notice that the transverse component of the current operatordoes not contribute enter into the continuity relation. We can therefore proceedas follows: (cid:68)(cid:104) ˆ n ( r , t ) , ˆ j (cid:96) ( r , t ) (cid:105)(cid:69) = − i (cid:90) β dλ (cid:68) ∇ · ˆ j ( r , j (cid:96) ( r , t − t + iλ ) (cid:69) . (C.7)Elementary vector analysis allows as to do the following manipulations: (cid:90) d r ∇ · ˆ j ( r , r , t ) = (cid:90) d r ∇ · (cid:104) ˆ j ( r , r , t ) (cid:105) − (cid:90) d r ˆ j ( r , · ∇ Φ( r , t ) . (C.8)In deriving the previous result, we have assumed that the external field vanishesat infinity and dropped the boundary term. We therefore reach (cid:90) d r ∇ · ˆ j ( r , r , t ) = (cid:90) d r j ( r , · E ( r , t ) . (C.9)Note that j ( r ,
0) is the longitudinal component of the current operator. Bykeeping in mind that ˆ j α ( r ,
0) is a longitudinal current component, we obtain J (1) (cid:96) ( r , t ) = (cid:88) α (cid:90) t −∞ dt (cid:90) d r (cid:90) β dλ (cid:68) ˆ j α ( r , j (cid:96) ( r , t − t + iλ ) (cid:69) E α ( r , t ) . (C.10)For two arbitrary ˆ A and ˆ B operators one can prove the following identity (cid:104) ˆ A (0) ˆ B ( t ) (cid:105) = 1 Z Tr (cid:104) e − β ˆ H ˆ A (0) ˆ B ( t ) (cid:105) = Tr (cid:104) ˆ ρe β ˆ H ˆ B ( t ) e − β ˆ H ˆ A (0) (cid:105) = (cid:104) ˆ B ( t − iβ ) ˆ A (0) (cid:105) . (C.11)We use the above identity in the following straightforward calculation i (cid:90) β dλ (cid:68) ˆ j α ( r , j (cid:96) ( r , t (cid:48) + iλ ) (cid:69) = (cid:90) t (cid:48) + iβt (cid:48) dτ (cid:68) ˆ j α ( r , j (cid:96) ( r , τ ) (cid:69) = (cid:90) ∞ t (cid:48) dτ (cid:68) ˆ j α ( r , j (cid:96) ( r , τ ) (cid:69) − (cid:90) ∞ t (cid:48) + iβ dτ (cid:68) ˆ j α ( r , j (cid:96) ( r , τ ) (cid:69) = (cid:90) ∞ t (cid:48) dτ (cid:68) ˆ j α ( r , j (cid:96) ( r , τ ) (cid:69) − (cid:90) ∞ t (cid:48) dτ (cid:68) ˆ j (cid:96) ( r , τ )ˆ j α ( r , (cid:69) = (cid:90) ∞ t (cid:48) dτ (cid:68)(cid:104) ˆ j α ( r , , ˆ j (cid:96) ( r , τ ) (cid:105)(cid:69) . (C.12)Eventually, for t (cid:48) = t − t , we can reach the following relation for the first-ordercurrent response J (1) (cid:96) ( r , t ) = − i (cid:88) α (cid:90) t −∞ dt (cid:90) d r (cid:90) ∞ t − t dτ (cid:68)(cid:104) ˆ j α ( r , , ˆ j (cid:96) ( r , τ ) (cid:105)(cid:69) E α ( r , t ) . (C.13)38ntroducing τ = τ and τ = t − t , we can rewrite the previous equation as J (1) (cid:96) ( r , t ) = − i (cid:88) α (cid:90) ∞ dτ (cid:90) ∞ τ dτ (cid:90) d r (cid:68)(cid:104) ˆ j α ( r , , ˆ j (cid:96) ( r , τ ) (cid:105)(cid:69) E α ( r , t − τ ) . (C.14)Using the following Fourier transformation E α ( r , t − τ ) = (cid:90) dωE α ( r , ω ) e − iω ( t − τ ) , (C.15)we arrive at J (1) (cid:96) ( r , ω ) = − i (cid:88) α (cid:90) ∞ dτ e iωτ (cid:90) ∞ τ dτ (cid:90) d r (cid:68)(cid:104) ˆ j α ( r , , ˆ j (cid:96) ( r , τ ) (cid:105)(cid:69) E α ( r , ω ) . (C.16)After performing an integration by parts, we arrive at J (1) (cid:96) ( r , ω ) = − i (cid:88) α (cid:40) e iωτ iω (cid:90) ∞ τ dτ (cid:90) d r (cid:68)(cid:104) ˆ j α ( r , , ˆ j (cid:96) ( r , τ ) (cid:105)(cid:69) (cid:12)(cid:12)(cid:12) τ →∞ τ → + (cid:90) ∞ dτ e iωτ iω (cid:90) d r (cid:68)(cid:104) ˆ j α ( r , , ˆ j (cid:96) ( r , τ ) (cid:105)(cid:69) (cid:41) E α ( r , ω )= (cid:88) α (cid:40) i (cid:90) ∞ dτ (cid:90) d r (cid:68)(cid:104) ˆ j α ( r , , ˆ j (cid:96) ( r , τ ) (cid:105)(cid:69) − i (cid:90) ∞ dτ e iωτ (cid:90) d r (cid:68)(cid:104) ˆ j α ( r , , ˆ j (cid:96) ( r , τ ) (cid:105)(cid:69) (cid:41) E α ( r , ω ) iω . (C.17)Using the linear-response definition given in Eq. (12), we find the followingrelation for the first-order current in the scalar potential gauge: J (1) (cid:96) ( r , ω ) = − (cid:88) α (cid:90) d r [ χ j (cid:96) j α ( r , r , ω ) − χ j (cid:96) j α ( r , r , E α ( r , ω ) iω . (C.18)Equivalently, in the wave-vector space, we have: J (1) (cid:96) ( q , ω ) = − (cid:88) q (cid:48) (cid:88) α [ χ j (cid:96) j α ( q , q (cid:48) , ω ) − χ j (cid:96) j α ( q , q (cid:48) , E α ( q (cid:48) , ω ) iω . (C.19)Comparing Eq. (63) with Eq. (C.19), we conclude that the following gauge-invariance identity must hold true: (cid:104) ˆ κ (cid:96)α ( q , q (cid:48) ) (cid:105) = lim ω → χ j (cid:96) j α ( q , q (cid:48) , ω ) . (C.20)For a longitudinal external field E α ( q (cid:48) , ω ) || q (cid:48) , the above relation is valid at fi-nite q , q (cid:48) . However, at finite q (cid:48) the response to a transverse field E α ( q (cid:48) , ω ) ⊥ q (cid:48)
39s not obtainable in the scalar potential gauge. The magnetic dipole coupling iscaptured only in calculations carried out in the vector potential gauge. There-fore, in the case of a transverse electric field, the above sum rule is not valid forfinite q (cid:48) . We note, however, that in the local q (cid:48) = 0 limit, results in the scalarand vector potential gauges are identical. Appendix D. Continuity relations for multi-photon current opera-tors
The field-dependent current operator is given byˆ J α ( r , t ; A ) = ˆ j α ( r , t ) + (cid:88) β (cid:90) d r (cid:48) ˆ κ αβ ( r , r (cid:48) , t ) A β ( r (cid:48) , t )+ (cid:88) βγ (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) ˆ ξ αβγ ( r , r (cid:48) , r (cid:48)(cid:48) , t ) A β ( r (cid:48) , t ) A γ ( r (cid:48)(cid:48) , t ) + . . . . (D.1)The charge conservation law (i.e. the continuity constitution law), i.e. d ˆ n/dt = i [ ˆ H + ˆ V , ˆ n ] = − ∇ · ˆ J , can be written more explicitly as following: − i ˆ H − (cid:88) β (cid:90) (cid:18)(cid:90) dλ ˆ J β ( r , t ; λ A ) (cid:19) A β ( r , t ) d r , ˆ n ( r , t ) = − ∇ · ˆ J ( r , t ; A ) . (D.2)We now expand the left-hand side of the above relation in terms of the vectorpotential (up to quadratic order):1 i (cid:104) ˆ H , ˆ n ( r , t ) (cid:105) − i (cid:88) β (cid:90) (cid:18)(cid:90) dλ ˆ J β ( r , t ; λ A ) (cid:19) A β ( r , t ) d r , ˆ n ( r , t ) = 1 i (cid:104) ˆ H , ˆ n ( r , t ) (cid:105) − i (cid:88) β (cid:90) d r (cid:48) (cid:104) ˆ j β ( r (cid:48) , t ) , ˆ n ( r , t ) (cid:105) A β ( r (cid:48) , t ) − i (cid:88) βγ (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) [ˆ κ βγ ( r (cid:48) , r (cid:48)(cid:48) , t ) , ˆ n ( r , t )] A β ( r (cid:48) , t ) A γ ( r (cid:48)(cid:48) , t ) + O ( A ) . (D.3)Similarly, we can expand the right-hand side of Eq. (D.2) as following (up toquadratic order): ∇ · ˆ J ( r , t ; A ) = ∇ · ˆ j ( r , t ) + (cid:88) αβ (cid:90) d r (cid:48) ∂∂r α ˆ κ αβ ( r , r (cid:48) , t ) A β ( r (cid:48) , t )+ (cid:88) αβγ (cid:90) d r (cid:48) (cid:90) d r (cid:48)(cid:48) ∂∂r α ˆ ξ αβγ ( r , r (cid:48) , r (cid:48)(cid:48) , t ) A β ( r (cid:48) , t ) A γ ( r (cid:48)(cid:48) , t ) + O ( A ) . (D.4)40he two relations given in Eq. (D.3) and Eq. (D.4) must be equal for an arbitraryvector potential. Equating the zero-th order term in the two expressions weobtain the one-photon continuity relation:1 i (cid:104) ˆ H , ˆ n ( r , t ) (cid:105) = (cid:88) α ∂∂r α ˆ j α ( r , t ) . (D.5)Similarly, equating the first- and second-order order terms, we obtain the two-and three-photon continuity relations, − i (cid:104) ˆ j β ( r (cid:48) , t ) , ˆ n ( r , t ) (cid:105) = (cid:88) α ∂∂r α ˆ κ αβ ( r , r (cid:48) , t ) (D.6)and − i [ˆ κ βγ ( r (cid:48) , r (cid:48)(cid:48) , t ) , ˆ n ( r , t )] = (cid:88) α ∂∂r α ˆ ξ αβγ ( r , r (cid:48) , r (cid:48)(cid:48) , t ) . (D.7)In the Schr¨odinger picture of time evolution, we have (cid:104) ˆ H , ˆ n ( r ) (cid:105) = − (cid:88) α (cid:18) − i ∂∂r α (cid:19) ˆ j α ( r ) (D.8) (cid:104) ˆ j β ( r (cid:48) ) , ˆ n ( r ) (cid:105) = (cid:88) α (cid:18) − i ∂∂r α (cid:19) ˆ κ αβ ( r , r (cid:48) ) , (D.9)and 12 [ˆ κ βγ ( r (cid:48) , r (cid:48)(cid:48) ) , ˆ n ( r )] = (cid:88) α (cid:18) − i ∂∂r α (cid:19) ˆ ξ αβγ ( r , r (cid:48) , r (cid:48)(cid:48) ) . (D.10)By plugging Eq. (D.9) in Eq. (D.10) we find12 (cid:104) [ˆ j γ ( r (cid:48)(cid:48) ) , ˆ n ( r (cid:48) )] , ˆ n ( r ) (cid:105) = (cid:88) αβ (cid:32) − i ∂∂r (cid:48) β (cid:33) (cid:18) − i ∂∂r α (cid:19) ˆ ξ αβγ ( r , r (cid:48) , r (cid:48)(cid:48) ) . (D.11)Let us just focus on the two-photon (diamagnetic) coupling which, in Fourierspace, reads [ˆ j β ( − q (cid:48) ) , ˆ n ( q )] = (cid:88) α q α ˆ κ αβ ( q , q (cid:48) ) . (D.12)Note that ˆ n ( r ) = (cid:88) q ˆ n ( q ) e i q · r , (D.13)ˆ j β ( r (cid:48) ) = (cid:88) q (cid:48) ˆ j β ( − q (cid:48) ) e − i q (cid:48) · r , (D.14)ˆ κ αβ ( r , r (cid:48) ) = (cid:88) q , q (cid:48) ˆ κ αβ ( q , q (cid:48) ) e i q · r e − i q (cid:48) · r (cid:48) . (D.15)41he Fourier transform of the density operator in first quantization readsˆ n ( q ) = (cid:88) i e − i q · ˆ r i , (D.16)where (cid:80) i runs over all particles. Similarly, the one-photon current operator isgiven by ˆ j β ( q ) = 12 (cid:88) i (cid:8) v β ( ˆ p i ) e − i q · ˆ r i + e − i q · ˆ r i v β ( ˆ p i ) (cid:9) , (D.17)where v β ( p ) = ∂ p β H ( p ) is the β -th component of the velocity operator and ˆ p i = − i ∇ i is the momentum operator corresponding to the i -th particle. Commutingcurrent and density operators as given in Eqs. (D.17) and (D.16) we then find:[ˆ j β ( − q (cid:48) ) , ˆ n ( q )] = (cid:88) ij [ v β ( ˆ p i ) , e − i q · ˆ r j ] e i q (cid:48) · ˆ r i + e i q (cid:48) · ˆ r i [ v β ( ˆ p i ) , e − i q · ˆ r j ]2= v β ( − q )ˆ n ( q − q (cid:48) ) . (D.18)Therefore, the particle-number conservation law leads to the following rela-tion v β ( − q )ˆ n ( q − q (cid:48) ) = (cid:88) α q α ˆ κ αβ ( q , q (cid:48) ) . (D.19)For the case of a homogenous electron liquid, by considering (cid:104) ˆ n ( q − q (cid:48) ) (cid:105) = nδ q , q (cid:48) with n the total particle density, we find nv β ( − q ) = (cid:88) α q α (cid:104) ˆ κ αβ ( q , q ) (cid:105) . (D.20)For a parabolic-band model, H ( p ) = p / m , we obtain: − q β nm = (cid:88) α q α (cid:104) ˆ κ αβ ( q , q ) (cid:105) , (D.21)which results in the well-known diamagnetic contribution (cid:104) ˆ κ αβ ( q , q ) (cid:105) = − nm δ αβ . (D.22)Note that in the main text we convert particle density/current to charge den-sity/current by including the electron charge: ˆ n → − e ˆ n , ˆ j α → − e ˆ j α , andˆ κ αβ → ( − e ) ˆ κ αβ . 42 ppendix E. Proof for (cid:80) q K = 0 In order to prove the constraint (cid:80) q K = 0, we first simplify the expressionfor K . From (75), we have K = (cid:88) α α [Π n ; α α ( Q , Ω ) + Π n ; α α ( Q , Ω )] [ q ,α A α ( q , ω )] − (cid:88) α [Π n ; nα ( Q , Ω ) ω A α ( q , ω ) + Π n ; nα ( Q , Ω )Φ( q , ω ) q ,α ]+ [Π n ; nn ( Q , Ω ) + Π n ; nn ( Q , Ω )] [ ω Φ( q , ω )] . (E.1)The definitions of Π n ; α α and Π n ; nn contain the intrinsic permutation symme-try operation. We can therefore simplify the above relation as follows K = 2 (cid:40) (cid:88) α α Π n ; α α ( Q , Ω ) [ q ,α A α ( q , ω )] + Π n ; nn ( Q , Ω ) [ ω Φ( q , ω )] (cid:41) − (cid:88) α [Π n ; nα ( Q , Ω ) ω A α ( q , ω ) + Π n ; nα ( Q , Ω )Φ( q , ω ) q ,α ] . (E.2)By plugging Eq. (76) in Eq. (E.2), we arrive at K = 2 (cid:88) α α Π n ; α α ( Q , Ω ) q ,α iω E α ( q , ω ) + 2 (cid:88) α Π n ; nα ( Q , Ω ) ω q ,α ω Φ( q , ω ) − (cid:88) α [Π n ; nα ( Q , Ω ) ω A α ( q , ω ) + Π n ; nα ( Q , Ω )Φ( q , ω ) q ,α ] . (E.3)We can collect the second and third terms in the above relation and reach K = (cid:88) α (cid:40) (cid:88) α Π n ; α α ( Q , Ω ) q ,α − Π n ; nα ( Q , Ω ) ω (cid:41) E α ( q , ω )+ i (cid:40) ω (cid:88) α Π n ; nα ( Q , Ω ) q ,α − ω (cid:88) α Π n ; nα ( Q , Ω ) q ,α (cid:41) Φ( q , ω ) . (E.4)By considering the definition of U α ( Q , Ω ), we can simplify Eq. (E.4) asfollowing: K = (cid:88) α U α ( Q , Ω ) E α ( q , ω ) + i (cid:88) α (cid:104) Π n ; nα ( Q , Ω ) q ,α ω − Π n ; nα ( Q , Ω ) q ,α ω (cid:105) Φ( q , ω ) ω ω . (E.5)43y using Eq. (79), we reach K = (cid:88) α (cid:20) Π n ; nα ( Q , Ω ) q ,α ω − Π n ; nα ( Q , Ω ) q ,α ω (cid:21) Φ( q , ω ) ω ω . (E.6)By considering Eq. (81) we can prove the following relation (cid:88) α Π n ; nα ( Q , Ω ) q ,α ω = 2 (cid:88) α α Π n ; α α ( Q , Ω ) q ,α q ,α ω ω = 2 (cid:88) α α Π n ; α α ( Q , Ω ) q ,α q ,α ω ω [Note: Relabelling: α ↔ α ]= 2 (cid:88) α α Π n ; α α ( Q , Ω ) q ,α q ,α ω ω [Note: Permutation symmetry]= (cid:88) α Π n ; nα ( Q , Ω ) q ,α ω . (E.7)By plugging the above relation in Eq. (E.6) we get K = 0 and therefore (cid:80) q K = 0. Following a very similar approach, one can prove that (cid:80) q L1