Gauge invariance and Ward identities in nonlinear response theory

We present a formal analysis of nonlinear response functions in terms of correlation functions in real- and imaginary-time domains. In particular, we show that causal nonlinear response functions, expressed in terms of nested commutators in real time, can be obtained from the analytic continuation of time-ordered response functions, which are more easily amenable to diagrammatic calculation. This generalizes the well-known result of linear response theory. We then use gauge invariance arguments to derive exact relations between second-order response functions in density and current channels. These identities, which are non-perturbative in the strength of inter-particle interactions, allow us to establish exact connections between nonlinear optics calculations done in different electromagnetic gauges.


Introduction
Nonlinear optical phenomena [1][2][3][4], arise from the response of an electronic system to intense electromagnetic fields. Going beyond linear response, many interesting optical effects can occur including harmonic generation, wave mixing, and saturable absorption [3]. As in the case of linear response theory, calculations of nonlinear response functions to an electromagnetic field can be carried 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) and the vector potential gauge (Φ = 0, A = 0). Nonlinear response functions are usually expressed in terms of complicated and lengthy relations, which are cumbersome to evaluate [2][3][4]. Therefore, one often needs to introduce some approximations in order to handle technical difficulties. For example, in the so-called "optical" limit, the photon momentum is neglected and one needs to study only the local response to a homogenous time-dependent electric field. This treatment is often called electric-dipole approximation [2,3]. In this approximation, the field-particle interaction can be written in two major gauges: H Φ = H(p) − eE(t) · r or H A = H(p + eA(t)), where H is the Hamiltonian in the absence of radiation, and r and p ≡ −i ∇ r stand for the position and momentum operators, respectively. The equivalence of these two gauges for an arbitrary physical observable in the electric-dipole approximation is well discussed in the literature [5][6][7][8][9].
For optical properties of insulators, we mainly need to consider inter-band transitions while the intra-band contribution is more relevant in metals and semimetals. To calculate the inter-band optical response, it is convenient to utilize the vector-potential gauge [10] in which one can simply drop the photon momentum from the very beginning of the calculation. Apart from pure intraand inter-band transitions, an extra class of transitions emerges only in the analysis of nonlinear response functions containing n-point correlation functions with n ≥ 3. This new contribution originates from mixed transitions which contain both intra-and inter-band processes [11,12]. Employing a generic density gauge, Φ(r, t), is more suitable for evaluating the contributions due to all intra-band, inter-band and mixed transitions in the optical limit [11][12][13]. This is because only a single non-interacting Feynman diagram is required to be evaluated in the density gauge at any order of perturbation theory, while in the vector potential gauge there are two, four, and eight diagrams that need to 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 beyond the electric-dipole approximation. If one keeps the photon momentum up to linear order in the nonlinear conductivity, it is possible to capture nonlocal effects in the electric-quadrupole and magnetic-dipole approximations [13,15]. These new contributions are very important for the second-order nonlinear response of inversion symmetric systems [3,16,17]. Note that in the scalar potential gauge, one loses the transverse contribution to the current, which originates from the magnetic-dipole moment [15]. In fact, as discussed in this Article, in order to calculate the nonlocal current it is preferable to perform the calculation in the vector 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 invariance in nonlinear response theory is devoted to exploring the gauge choice (the two specific aforementioned gauges in the electric-dipole approximation) in noninteracting atoms [18] and crystalline solids [8]. To the best of our knowledge, a generic analysis of gauge invariance for a spatially inhomogeneous external field in nonlinear response theory as applied to interacting electron systems is still unexplored and it is one of the main motivations of this work. Our analysis of nonlinear response functions covers all local (electric-dipole) and nonlocal (other electric and magnetic-multipoles) effects and can be seen as a generalization of the gauge invariance analysis in the context of linear response theory, which is textbook material [19][20][21]. Nonlocal response functions are also of interest when the momentum exchanged between photons (external field) and charged particles cannot be ignored. This is the case of the photon-drag effect [22][23][24] and of other effects where the photon momentum needs to be taken into account. For example, our analysis may be relevant to analyze experimental results of spectroscopy based on free-electron lasers [25,26].
Gauge invariance, i.e. the independence of physical observables on the electromagnetic gauge one chooses for calculations, imposes severe constraints on the theory of linear response [19][20][21]. Indeed, 1) the fact that a static and purely longitudinal vector potential cannot produce any physical current implies the vanishing of the longitudinal current response at any finite wave vector q; 2) similarly, the fact that a static and quasi-homogeneous vector potential cannot produce any physical current implies the vanishing of the longitudinal and transverse current-current response functions for q → 0. More generally, gauge invariance imposes precise relationships between linear response functions to scalar and vector potentials.
In this Article, we first present a formal analysis of the nonlinear response theory in the real and imaginary time domains. We show that causal nonlinear response functions, expressed in terms of nested commutators in real time, can be obtained from the analytic continuation of time-ordered response functions, which are more easily amenable to diagrammatic calculation. This generalizes the well-known result of linear response theory.
We then provide a theoretical study of the impact of gauge invariance on the second-order response functions. We finally report Ward identities which must be fulfilled in order to ensure gauge invariance at any order of perturbation theory. In Section 2, we report an explicit definition of nonlinear response functions. In Section 3, we discuss the spectral representation of second-order response functions in real-and imaginary-time domains. In Section 4, we consider a specific kind of external field, i.e. electromagnetic radiation treated classically. In Sections 5 and 6, we explicitly discuss gauge invariance for linear and secondorder response functions, respectively. Finally, in Section 7, we present a set of Ward identities for nonlinear response functions.

Nonlinear response theory in real time
We consider a many-body system in thermal equilibrium which is described by a Hamiltonian denoted byĤ. We then turn on an external field, which can be modelled by a field-particle interaction term denoted byV . We use theŜmatrix approach [27] to consider the effect of the interaction in a perturbative way. In this approach, one has an adiabatic time evolution of the wave function from an unperturbed state at time t = −∞, i.e. |ψ 0 , to the perturbed state at time t. Therefore, we have |ψ(t) =Ŝ(t, −∞)|ψ 0 where theŜ-matrix is given by (setting = 1 for the sake of simplicity) Here, T stands for the time-ordering operation andV (t) is the perturbative part of the Hamiltonian in the interaction representation, i.e.V (t) = e itĤV e −itĤ . In a similar way, any observable evolves in the interaction picture aŝ Here,Â denotes the second-quantized representation of the A operator. The expectation value ofÂ at time t is given by Thanks to the Dyson expansion, we have (see Appendix A) where Θ(τ ) is the Heaviside function. The field-particle interaction can be formally written as followŝ whereB i (t) indicates linear coupling whileĈ ij (t) andD ijk (t) stand for the nonlinear couplings. We can visualize these linear and nonlinear coupling terms as vertices in the diagrams shown in Fig. 1. For the case of light-matter interactions, we term the B, C, and D vertices as single-photon, two-photon, and three-phonon vertices, respectively. Here, we first assume a linear coupling viâ B i to the external fields F i and later generalize our formal theory to the nonlinear couplingsĈ ij (t) andD ijk (t). Note that vertex couplings likeĈ ij (t) and D ijk (t) only emerge in the non-relativistic limit of a low-energy theory because of the nonlinear energy dispersion of the band Hamiltonian ( k ∼ n≥0 n k n ) and they do not exist in quantum electrodynamics [28]. After neglecting non-B iĈ ijDijk Figure 1: Diagrammatic representation of the linear and nonlinear couplings as vertices (solid dots) 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 The previous expression leads us to define the linear retarded response functions [21] χ ABi (t; the second-order retarded response functions and the third-order ones Of course, one can simply extend these definitions to higher-order response functions. The symbol P in Eqs. (8) and (9) is there to ensure a permutation symmetry. Since the quantity Π m j=1 F j (t−τ j ) in Eq. (6) is symmetric with respect to 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 response function. For instance, following Ref. [3], one can decompose the m-th order response function in the sum of symmetric and anti-symmetric contributions, i.e. χ = χ s + χ a with and It is obvious that χ s ( χ a ) is symmetric (anti-symmetric) with respect to the permutation between (B 1 , τ 1 ) and (B 2 , τ 2 ). Since the term Π m j=1 F j (t − τ j ) is symmetric for all orders of permutation, the anti-symmetric part plays no role in determining the expectation value ofÂ. The only non-vanishing contribution to the latter originates from completely symmetric response functions. This is the reason why in the definition of the m-th order nonlinear response function we have introduced the sum P = 1 m! P , where P stands for the sum over all permutations among the dummy variables (B i , τ i ). This kind of symmetry is usually called intrinsic permutation symmetry [3].
By using the cyclic properties of the trace, we can then eliminate "t" in Eqs. (7), (8), and (9): The invariance under time translations of the previous response functions is evident. They just depend on the time differences τ i = t − t i . This can be powerfully used by Fourier transforming to the frequency domain: Here, ω Σ = i ω i . As usual [21], we have assumed that all frequencies contain an infinitesimal positive imaginary part η i → 0 + . This stems from the need to fit periodic perturbations F i (t) = F i (ω)e −iωt e −ηt + c.c. into the response theory formalism, making sure that these vanish in the remote past. In the frequency 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. Ô = Tr[ρÔ] = λ P λ λ|Ô|λ whereρ is the density matrix, P λ = λ|ρ|λ = e −βE λ /Z, Z = λ e −βE λ is the partition function, and β = 1/(k B T ) with T the temperature. Notice that |λ and E λ are the exact eigenstates and eigenvalues of the many-body HamiltonianĤ in the absence of radiation. The same 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 this nonlinear case too.
With the definitions given in Eqs. (12), (13), and (14), we can take advantage of diagrammatic techniques to evaluate nonlinear response functions. In Fig. 2, we give a diagrammatic representation of the m-th order response function in both time and frequency domains. The generalization of the linear-coupling theory to the nonlinear couplings is straightforward as it is simply based on the inclusion of higher-order photon vertices in the Feynman diagrams.

Spectral representation and analytic continuation
In an interacting many-body system, it is not easy to handle all the diagrams for linear and nonlinear response functions. It is often convenient to do the calculations with time-ordered response functions in imaginary time and then come back to real time by an analytic continuation [27]. The imaginary-time representation of response functions, also known as "Matsubara representation", is also extremely useful for the finite temperature analysis. In order to derive the connection between causal (real-time) response function and time-ordered response functions in imaginary time, we proceed in three steps. First we derive the so-called spectral representation (also known as "Lehmann representation" or "exact eigenstates representation") for second-order response functions. Next, making use of this representation, it is shown that the second-order causal response functions in real time are connected to time-ordered second-order response functions in imaginary time by the well-known analytic continuation procedure [27]. Finally, reasoning by induction, we prove that the analytic continuation procedure works for response functions of arbitrary order.

Real time (frequency) spectral representation: second-order response
According to Eq. (13), the second-order retarded response function in the real-time representation reads as follows: where Performing the thermodynamic average, we obtain We exchange λ 1 , λ 2 and λ 3 dummy labels in the second term of the above relation in order to the collect it with the first term. We use a similar trick and collect the last two terms together. Consequently, we arrive at Notice that λ 3 |Â(0)|λ 1 = A λ3λ1 and P λ1λ2 = P λ1 − P λ2 . We extract the time-dependent part as follows where E λ1λ2 = E λ1 − E λ2 . Using Eq. (16), we obtain the frequency-domain representation of the second-order response function: Considering Eq. (22), we can easily perform the integrations in Eq. (23) and we reach Eventually, by considering the intrinsic permutation, we arrive at the following final expression for the retarded second-order response function: We should set η 1 = η 2 in order to turn on all the external fields at the same rate in the adiabatic switching-on procedure.

Imaginary time (frequency) spectral representation: second-order response
In the imaginary time, τ ≡ it, the time evolution of an operator obeys: O(τ ) = e τ HÔ e −τ H . We should notice that only in this subsection τ represents an imaginary time. We generalize the periodic property of a bosonic correlation function, i.e. f (τ ) = f (τ + β), [27] to a multi-variable correlation function as f (τ 1 , τ 2 , . . . , τ m ). In other words, for an arbitrary choice of τ i , we have f (τ 1 , τ 2 , . . . , τ i , . . . , τ m ) = f (τ 1 , τ 2 , . . . , τ i + β, . . . , τ m ) . Therefore, we can introduce the following Fourier and inverse-Fourier transformations: where ν i = 0, ±2π/β, ±4π/β . . . stands for the bosonic Matsubara frequency (energy). Now, we employ the Matsubara representation for the case of the second-order response in order to generalize the first-order Matsubara technique to higher-order response functions. The imaginary time-ordered second-order response function is defined by After performing the time-ordering operation, we arrive at Since 0 ≤ τ i ≤ β-see Eq. (26)-we have Θ(−τ i ) = 0 and Θ(τ i ) = 1 which implies The above equation can be rewritten in a compact form as follows where P is the usual intrinsic permutation operation symbol. The thermodynamic average can be taken in the spectral representation: According to the time-evolution convention and considering Eq. (26) for the Fourier transformation definition in the imaginary-time domain, we Fourier transform χ AB1B2 (τ 1 , τ 2 ) obtaining Using that P λ = e −βE λ /Z, we can simplify the above relationship as follows We can write the fraction in the first line of the previous equation in the following way (34) Using the above identity, we can simplify the expression of χ AB1B2 (iν 1 , iν 2 ) even more and reach After exchanging λ 2 with λ 3 , we obtain (36) By considering the intrinsic permutation operation, we reach the following compact form for the spectral representation of the second-order response in the Matsubara frequency domain: (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, we should set η 1 = η 2 in order to turn on all the external fields at the same rate in the adiabatic switching-on procedure.

Analytic continuation for higher order nonlinear response functions
We define the m-th order time-ordered response function as follows: where τ i denotes an imaginary time. Notice that the (−1) m /m! factor is required in order to achieve consistency with the real-time picture, see e.g. Eq. (25). The intrinsic permutation symmetry is implicitly taken into account through the presence of the time-ordering operation, T , and the 1/m! pre-factor is there to avoid multiple counting. A similar relation for the imaginary-time nonlinear correlation function was first reported in Ref. [29]. The response function can be expressed in terms of a "universal" kernel X (n) in the following manner: The universal kernel, which entirely accounts for the frequency dependence of the response, can be evaluated in a recursive manner as described below. In the zero-th order it is given by the statistical occupation factor: The first-order case then follows Following the derivation given above regarding the analytical continuation of the second-order response function, we find Similarly, for the third order we obtain and thus eventually for the n-th order we have (44) This is the desired recursion relation. Having the universal X (n) response function, one can calculate the n-th order physical response function after properly incorporating the form-factor part shown in the square bracket in Eq. (39) and summing on all degrees of freedoms λ i . We have already established that, up to second-order, the analytic continuation from time-ordered to causal response is effected by the replacement, iν i → ω i + iη i . The recursion relation (44) shows that the same procedure will also work for the n-th order response if it works for the n − 1-th response. Thus, we have provided an inductive proof of the analytic continuation procedure at all orders. Our proof is considerably simpler than the one first reported in Ref. [29].

Light-matter interaction: gauge transformation
After the above formal analysis of nonlinear response functions, we now proceed to discuss the issue of gauge invariance in nonlinear response theory, analysing the role of an external electromagnetic field. We then generalize the linear Ward identity stemming from gauge invariance and charge conservation to higher-order response functions [15,28].
In the interaction representation, the light-matter interaction can be written in terms of vector, A(r, t), and scalar, Φ(r, t), potentials in the space-time domain, (r, t): whereĴ α1 (r, t; A) is the charge current operator, which, in general, may depend on the electromagnetic vector potential A. The integration over the real parameter λ in the first term of Eq. (45) guarantees the fundamental relation between the current and the light-matter interaction Hamiltonian: for any value of A. Note that the time dependence of the operatorsV andĴ α originates from two sources: the interaction-picture of the time evolution, see Eq. (2), and the explicit time dependence of the external potentials. The time dependence of the charge-density operatorn stems from the interaction-picture of the time evolution. The charge-current operator in the interaction picture is given bŷ whereĤ tot =Ĥ+V . Here, the one-photon current component in the interaction picture is given bŷ For the two-photon current coupling we havê Similarly, three-and four-photon current couplings read as following: It is common to dubĵ α1 "paramagnetic" current operator andκ α1α2 ,ξ α1α2α3 , andζ α1α2α3α4 "multi-photon" current operators (which can have either a diamagnetic or a paramagnetic contribution to total current depending on the ground-state phase and the order of perturbation; jargon stemming from the theory of superconductivity [19]). Note, finally, that the time dependence of j α1 (r, t) (and all multi-photon current operators) stems from the interactionpicture of the time evolution, see Eq.
(2). Formally, we have the following perturbative series for the macroscopic charge density and current We now apply a gauge transformation to obtain gauge-invariant relations for the first-and second-order density and current. The gauge transformation is defined by the following map where Λ(r, t) is an arbitrary smooth function of r and t. Similarly, in Fourier space: The gauge transformation preserves the electric and magnetic fields, where we have set c = 1 for the speed of light in vacuum. In Fourier space, the field components read where αβγ is the Levi-Civita symbol. From now on, we drop the imaginary part of ω in the notation but keep in mind that ω ≡ ω + i0 + .
According to gauge invariance, the difference between physical observables calculated in two arbitrary gauges must be zero. This implies:

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 charge density and current: Note that our system is not assumed to be translationally invariant and consequently q and q are two independent variables. According to gauge invariance and using Eq. (55), the difference between the physical charge densities and currents calculated in two different gauges must be identically zero: In order to fulfill these constraints, the following gauge-invariance identities must be satisfied After implementing the above identities in Eq. (58) and using Eq. (57), we can write down the following gauge-invariant form of the linear density and current: The surprising absence of the magnetic field in these equations is explained by noting that the finite frequency components of the magnetic field are connected to the electric field by the two "internal" Maxwell equations, e.g. Faraday's law q × E(q, ω) = ωB(q, ω) and q · B(q, ω) = 0, which follow from Eq. (57).
If we expand the conductivity up to linear order in q we can generate terms proportional to the magnetic field. Because of the above relation for the first-order current, we have the following expression for the linear conductivity: Utilizing the diamagnetic sum rule for a longitudinal external field E(q , ω)||q , we have (see Ref. [21] and also Appendix C) For longitudinalĵ α : For a transverse field E(q , ω) ⊥ q the total static current response is finite and therefore the diamagnetic sum rule is valid only at q = 0: The diamagnetic sum rule identity is explicitly proven in Appendix C by using the well-known Kubo identity [30,31] and the continuity relation. This cancellation clarifies the importance of the diamagnetic component of the current operator, i.e. ακ α A α , in the linear conductivity. Finally, we recall that in the homogeneous (translational invariant) electron liquid we have q = q . In the clean system, the Drude divergence of the dc limit, ω → 0, is achieved in the reverse order of limits that is first setting q → 0 and then ω → 0. In this case, the paramagnetic response function vanishes, lim ω→0 lim q→0 χ j jα → 0, in a perfectly homogeneous electron liquid owing to momentum conservation [21]. As a consequence, the diamagnetic contribution leads to the Drude divergence at zero frequency.
The continuity constitution laws for the two-and three-photon diamagnetic current couplings read (see Appendix D): and A similar strategy can be employed to resolve a connection between the equaltime nested commutation of the paramagnetic current and the density operator with higher-order n-photon diamagnetic current operators. Let us focus on the two-photon diamagnetic current coupling, which, in Fourier space, obeys It is straightforward to evaluate [ĵ β (−q ),n(q)] and obtain (see Appendix D) Notice that e > 0 is the fundamental electric charge and v β (q) = ∂ q β H(q) is the velocity operator.
There is a simpler approach to prove the above identity. Using the scalar potential 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 calculated in the two gauges must be equal owing to gauge invariance. Therefore, we obtain Accordingly, we find the same identity given in Eq. (82) after cancelling the scalar potential from the above relation. The above relation is an identity which holds true to ensure the gauge invariance of the second-order density N (2) (q, ω Σ ). By using Eqs. (76) and (81) and considering the intrinsic permutation symmetry of Π n;α1α2 , one can prove that q2 K 1 = 0 is also fulfilled. This proof is reported in Appendix E.
We now proceed to obtain a gauge-invariance relation for the second-order current. The second-order current in the frequency and wave-vector domain can be written as follows (see Appendix B) where Π ;nn (Q 12 , Ω 12 ) = P χ jαnn (Q 12 , Ω 12 ) and Π ;nα2 (Q 12 , Ω 12 ) = χ j njα 2 (Q 12 , Ω 12 ) + χ j jα 2 n (Q 21 , Ω 21 ) For the following two correlation functions, one needs to be careful in the symmetrization process, which should be carried out as follows After performing the gauge transformation, we reach the following relation for the second-order gauge-induced current change, δJ (2) (q, ω Σ ): where L 2 (Q 12 , Ω 12 ) = α1α2 Π ;α1α2 (Q 12 , Ω 12 )q 1,α1 q 2,α2 − α2 Π ;nα2 (Q 12 , Ω 12 )ω 1 q 2,α2 and Because of gauge invariance, both q2 L 1 and L 2 must be identically zero, a fact that leads to the following relation for L 2 = 0 : Replacing Eq. (92) in Eq. (84) we find where V α2 (Q 12 , Ω 12 ) = 2 α1 Π ;α1α2 (Q 12 , Ω 12 ) q 1,α1 − Π ;nα2 (Q 12 , Ω 12 ) ω 1 . Eq. (93) is not a fully gauge-invariant relation for the second-order current because it contains the scalar potential Φ. This is because we have not yet used the condition q2 L 1 = 0. However, for practical reason we instead perform another gauge transformation and later one can show that q2 L 1 = 0 is also satisfied. Accordingly, after performing a gauge transformation to the above relation we obtain Because of the gauge invariance property of the physical charge current, we obtain another gauge invariance identity: By considering this new identity, we reach the following gauge-invariant relation of the second-order current Therefore, the second-order conductivity reads as follows: Since Eq. (95) is valid for an arbitrary electric field, we conclude that V α2 (Q 12 , Ω 12 ) = 0. Using this relation, we can simplify Eq. (92) as follows The same result can be obtained in an intuitive way similar to Eq. (83). The above equation represents an important Ward identity, which ensures the gauge invariance of the second-order current. Similarly to the case of the second-order density, it is easy to prove that q2 L 1 = 0-this result can be obtained by performing a calculation similar to that reported in Appendix E for the proof of q2 K 1 = 0.
Finally, we visualize linear, second-order, and third-order response functions in the scalar and vector potential gauges in Figs. 3 and Fig. 4, respectively. As it 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 response functions, i.e. Π ;α1...αm , respectively.

Continuity relation and Ward identities
In the absence of external fields there is no current flow (unless the ground state carries it because of e.g. spontaneous breakdown of time-reversal symmetry). We therefore have ∂ n(r, t) /∂t = 0 and ĵ α (r, t) = 0. In this regard, the continuity equation can be written order-by-order in perturbation theory in the following manner: In the frequency domain we therefore have: Replacing the above relation in Eq. (61), we arrive at We therefore recover the well-known linear-response Ward identity, where In a translationally-invariant system q = q and χ J J (q, q , ω) represents the longitudinal current-current response function.
Similarly to the linear-response case, we can obtain a "second-order Ward identity" from the second-order continuity equation. Using Eqs. (80), (96) and Eq. (100), we arrive at q1q2 α1α2 q Π ;α1α2 (Q 12 , Ω 12 ) − ω Σ Π n;α1α2 (Q 12 , Ω 12 ) We therefore conclude that Π n;α1α2 (Q 12 , Ω 12 ) = q Π ;α1α2 (Q 12 , Ω 12 )/ω Σ . Substituting this relation in Eq. (82), we obtain the following second-order Ward identity: 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: Notice that in a translationally invariant system, we have q = q Σ = m i=1 q i . Moreover, the m-th order conductivity, the desirable gauge-invariant response function, reads σ (m) α1...αm (q, q 1 , . . . , q m , ω 1 , . . . , ω m ) = Π ;α1...αm (q, q 1 , . . . , q m , ω 1 , . . . , ω m ) It is useful to discuss the relation between nonlinear conductivities and nonlinear density response functions. Considering Eq. (107) and Eq. (108) we find that Π n;...n (q, q 1 , . . . , q m , ω 1 , . . . , This relation provides a gauge-covariant prescription in order to relate the dynamical non-local conductivity to the density response function at each order of perturbation theory. In the optical (electric-dipole) approximation we can neglect the wave-vector dependence of the conductivity and therefore the wave-vector expansion of the density response function in the electric-dipole approximation reads as following: In a translationally-invariant system we have q = q Σ . The latter is therefore not an independent variable, as indicated by the notation "q Σ, [q 1,α1 . . . q m,αm ]". The quantity σ We can go beyond the electric-dipole approximation and consider electricquadrupole and magnetic dipole contributions by retaining the wave-vector dependence of the conductivity up to linear order: The above relation was used [13] in studying nonlinear plasmonic effects in graphene, which is a centro-symmetric material. Because of the latter property, higher multipole contributions are important in this material [13].

Summary
In this Article, we have used the equilibriumŜ-matrix approach in order to set up a theory that allows the calculation of nonlinear response functions in both real-and imaginary-time domains. We have discussed an analytical continuation procedure for the nonlinear response functions, which provides a prescription in order to obtain the retarded response function from the Matsubara one.
A large fraction of this work was devoted to the analysis of gauge invariance in the context of the nonlinear response theory. An inhomogeneous vector potential gauge is a complete gauge to take into account contributions from all multipoles, such as electric-dipole, electric-quadruple, and magnetic-dipole contributions. As a result of gauge invariance, a set of nonlinear Ward identities are obtained.

Appendix A. Perturbation theory
The macroscopic value of a generic operatorÔ(t) reads as following: Thanks to the Dyson expansion, we havê where we define D m 1 as follows Using Eq. (A.1), the macroscopic value ofÔ(t) is given by In order to simplify the terms proportional to i 2 in the above relation, we need to rewrite D 1 1V (t 1 )Ô(t)D 1 1V (t 1 ) in a proper form. To do so, we divide the square domain of the integral into two triangular domains (see Fig. A.5). We therefore find Using the above identity, we reach the following compact form for the secondorder perturbative terms in Eq. (A.4): Similarly, one can find the following simplified expression for the third-order perturbation Eventually, after plugging Eqs. (A.6) and (A.7) in Eq. (A.4), we arrive at the following relation for O(t): Using the Heaviside function Θ(τ ), we reach

Appendix B. Position and real-time representation of the linear and second-order responses
First of all, we calculate the following terms which contain different perturbative contributions. The average of the current operator reads Ĵ (r, t; A) = ĵ (r, t) + dr α1 κ α1 (r, r , t) A α1 (r , t) + dr dr α1α2 ξ α1α2 (r, r, r , t) A α1 (r , t)A α2 (r , t) + . . . .
and similarly, ξ α1α2 (r, r , r , t) = ξ α1α2 (r, r , r ) . Eventually, we have The other required terms are the following commutators in whichV (t) is defined in Eq. (45) as the light-matter interaction. Moreover, we calculate the following necessary term In a similar way we can write down the corresponding expressions for V (t 1 ),n(r, t) and V (t 2 ), V (t 1 ),n(r, t) , which are necessary to obtain the first-and second-order densities.
In a similar way, the second-order current in position space and time domain can be written as follows (B.8) Using Eqs. (12) and (13), we rewrite the previous relation as follows (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 final expressions 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: (B.10) Using the linear-response definition given in Eq. (12), we rewrite the previous relation as follows The second-order density is given by By using the definitions given in Eqs. (12) and (13), we arrive at 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 expressions for 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, we perform a Fourier transformation. For a general expression of the type (B.14) the Fourier-transformed counterpart reads as following: We notice that for the case of translationally-invariant system we have q = q Σ = m j 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-order current: By using the cyclic property of the trace, we can prove that: n(r 1 , t 1 ),ĵ (r, t) = Tr ρ n(r 1 , t 1 ),ĵ (r, t) = Tr [ρ,n(r 1 , t 1 )]ĵ (r, t) Notice that in the last equality of the above relation, we have used the Kubo identity [30,31]: By performing a straightforward calculation one can show that n(r 1 , t 1 ),ĵ (r, t) = It is important to notice that the transverse component of the current operator does not contribute enter into the continuity relation. We can therefore proceed as follows: Elementary vector analysis allows as to do the following manipulations: In deriving the previous result, we have assumed that the external field vanishes at infinity and dropped the boundary term. We therefore reach Note that j(r 1 , 0) is the longitudinal component of the current operator. By keeping in mind thatĵ α (r 1 , 0) is a longitudinal current component, we obtain We use the above identity in the following straightforward calculation Eventually, for t = t − t 1 , we can reach the following relation for the first-order current response dτ ĵ α (r 1 , 0),ĵ (r, τ ) E α (r 1 , t 1 ) .

(C.13)
Introducing τ 2 = τ and τ 1 = t − t 1 , we can rewrite the previous equation as (C.14) Using the following Fourier transformation we arrive at (C. 16) After performing an integration by parts, we arrive at Using the linear-response definition given in Eq. (12), we find the following relation for the first-order current in the scalar potential gauge: J (1) (r, ω) = − α dr 1 [ χ j jα (r, r 1 , ω) − χ j jα (r, r 1 , 0)] E α (r 1 , ω) iω . (C. 18) Equivalently, in the wave-vector space, we have: Comparing Eq. (63) with Eq. (C. 19), we conclude that the following gaugeinvariance identity must hold true: κ α (q, q ) = lim ω→0 χ j jα (q, q , ω) . (C.20) For a longitudinal external field E α (q , ω)||q , the above relation is valid at finite q, q . However, at finite q the response to a transverse field E α (q , ω) ⊥ q is not obtainable in the scalar potential gauge. The magnetic dipole coupling is captured only in calculations carried out in the vector potential gauge. Therefore, in the case of a transverse electric field, the above sum rule is not valid for finite q . We note, however, that in the local q = 0 limit, results in the scalar and vector potential gauges are identical.

Appendix D. Continuity relations for multi-photon current operators
The field-dependent current operator is given bŷ where i runs over all particles. Similarly, the one-photon current operator is given byĵ where v β (p) = ∂ p β H(p) is the β-th component of the velocity operator andp i = −i∇ i is the momentum operator corresponding to the i-th particle. Commuting current and density operators as given in Eqs. For the case of a homogenous electron liquid, by considering n(q −q ) = nδ q,q with n the total particle density, we find Note that in the main text we convert particle density/current to charge density/current by including the electron charge:n → −en,ĵ α → −eĵ α , and κ αβ → (−e) 2κ αβ .