Spectral Density and Sum Rules for Second-Order Response Functions
Barry Bradlyn Peter Abbamonte







Matlab















Γ



α
Γ



    bbradlyn@illinois.eduDepartment of Physics, University of Illinois at Urbana-Champaign, Urbana IL 61801, USAAnthony J. 
    Leggett Institute for Condensed Matter Theory, University of Illinois at Urbana-Champaign, Urbana IL 61801, USADepartment of Physics, University of Illinois at Urbana-Champaign, Urbana IL 61801, USAMaterials Research Laboratory, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
    Sum rules for linear response functions give powerful and experimentally-relevant relations between frequency moments of response functions and ground state properties. 
    In particular, renewed interest has been drawn to optical conductivity and density-density sum rules and their connection to quantum geometry in topological materials. 
    At the same time, recent work has also illustrated the connection between quantum geometry and second-order nonlinear response functions in quantum materials, motivating the search for exact sum rules for second-order response that can provide experimental probes and theoretical constraints for geometry and topology in these systems. 
    Here we begin to address these questions by developing a general formalism for deriving sum rules for second-order response functions. 
    Using generalized Kramers-Kronig relations, we show that the second-order Kubo formula can be expressed in terms of a spectral density that is a sum of Dirac delta functions in frequency. 
    We show that moments of the spectral density can be expressed in terms of averages of equal-time commutators, yielding a family of generalized sum rules; furthermore, these sum rules constrain the large-frequency asymptotic behavior of the second harmonic generation rate. 
    We apply our formalism to study generalized f-sum rules for the second-order density-density response function and the longitudinal nonlinear conductivity. 
    We show that for noninteracting electrons in solids, the generalized f-sum rule can be written entirely in terms of matrix elements of the Bloch Hamiltonian.  
    Finally, we derive a family of sum rules for rectification response, determining the large-frequency asymptotic behavior of the time-independent response to a harmonic perturbation.
    Spectral Density and Sum Rules for Second-Order Response Functions
        Peter Abbamonte
        April 29, 2024
    ==================================================================
    
    
    
    § INTRODUCTION
    
    Sum rules for response functions represent some of the few exact results in many-body physics <cit.>.
    They are of profound fundamental importance, and are critical for interpreting spectroscopy experiments where they are used to calibrate the absolute scale of measurements <cit.>.
    The best known sum rule, known as the f-sum rule, relates the total electron density in a solid to the absorption probability of light, placing rigorous constraints on optical properties of materials.
    Furthermore, the f-sum rule is of practical utility to experimentalists when it comes to extrapolating the frequency dependence of measured response functions <cit.>.
    The related conductivity sum rule analogously determines the frequency of plasma oscillations in terms of the frequency-integrated AC susceptibility of a solid.
    Via Kramers-Kronig relations, sum rules also tell us about the large-frequency behavior of response functions <cit.>.
    
    Recently, attention has been paid to sum rules derived for more general moments of the conductivity tensor. 
    In particular, Ref. <cit.> derived a sum rule relating the first negative moment of the optical (AC) conductivity tensor to polarization fluctuations in the ground state of a many body system. 
    For noninteracting insulators, the polarization fluctuations can be expressed in terms of the quantum geometry of the occupied electronic states via the Fubini-Study metric, which has proved important in determining the superfluid weight in flat band systems <cit.>. 
    More generally, for interacting systems the polarization fluctuation sum rule can be expressed in terms of the “quantum weight” that generalizes the Fubini-Study metric to interacting systems <cit.>. 
    The connection between sum rules and quantum geometry has brought renewed interest to the study of linear response functions.
    
    However, outside of specialized applications to nonlinear optics <cit.>, the main focus of study has been sum rules for linear response functions.
    Notable exceptions are Refs. <cit.> which considered a family of nonlinear conductivity sum rules derived from the instantaneous current response to an applied (spatially uniform) electric field, as well as Ref. <cit.> which derives a sum rule for the nonlinear Hall conductivity. 
    However, the resurgence of interest in nonlinear—and in particular second-order—response functions in the study of topological materials <cit.> motivates the study of nonlinear sum rules more systematically. 
    For instance, Ref. <cit.> has shown that the nonlinear conductivity of quantum systems is connected to the Riemannian geometry of quantum states. 
    Sum rules for nonlinear response functions could thus place important constraints on, and provide experimental probes for, quantum geometry in insulators. 
    They would also enable experimenters to properly calibrate nonlinear spectroscopy measurements in the same manner customary for linear spectroscopies.
    
    In this work, we will take the first steps towards addressing these questions by studying sum rules for second-order response functions. 
    We will ultimately derive a set of sum rules for second harmonic generation response of a system to general perturbation. 
    First, in Sec. <ref> we will review the Kubo formula for second-order response functions. 
    We will introduce a representation for the response function in terms of a spectral density, generalizing the Kramers-Kronig relations of linear response theory. 
    Next, in Sec. <ref> we will show that frequency moments of the spectral density determine the large-frequency behavior of the second-order response function. 
    Furthermore, we will show that these frequency moments can be expressed in terms of averages of equal-time commutators, thus furnishing a set of second-order sum rules. 
    We will show that by focusing on second harmonic generation response, we arrive at a sum rule relevant to experiments. 
    In Sec. <ref> we will specialize to consider the second-order density response function χ_ρρ^(2), which quantifies the response of the local charge density ρ(𝐫) to a scalar potential perturbation at second order. 
    Applying our sum rules to χ_ρρ^(2) allows us to derive a generalization of the f-sum rule relevant to nonlinear spectroscopy. 
    This sum rule could enable determination of absolute cross sections for nonlinear optical processes such as second harmonic generation, parametric down conversion, or optical rectification.
    As an example of our formalism, in Sec. <ref> we evaluate our nonlinear f-sum rule for a system of free fermions in a periodic potential. 
    By exploiting charge conservation, we show in Sec. <ref> that our nonlinear f-sum rule also implies a sum rule for the longitudinal component of the nonlinear conductivity tensor, applicable to optical experiments. 
    Finally, in Sec. <ref> we turn our attention to second-order rectification, the time-independent response to a harmonic perturbation. 
    We show that for operators that can be written as time derivatives (such as the uniform current, which is the time derivative of the position operator), the large frequency asymptotic behavior of the rectification response is determined by a family of sum rules 
    We conclude in Sec. <ref> with a discussion of potential experimental and theoretical applications. 
    
    For notational convenience, we will work in units where ħ=c=e=1. 
    Additionally, for integrals over the domain (-∞,∞) we will leave the limits of integration implicit.
    
    
    
    § SECOND-ORDER RESPONSE AND THE SPECTRAL DENSITY
    
    
    Our setting for studying nonlinear response is a system with unperturbed Hamiltonian H_0, coupled to an external field. 
    We model the coupling to the external field by the perturbing Hamiltonian
    
        H_1(t)=e^ϵ tf(t)B,
     
    where f(t) is the amplitude of the externally applied field, B is the operator to which the field couples, and ϵ is a small positive infinitesimal that we will take to zero at the end of any calculation (it ensures that the perturbation goes to zero infinitely far in the past). 
    We will be interested in how the average ⟨ A⟩(t) of some observable A depends on the perturbation f(t). 
    We can expand in powers of f(t) to write
    
        ⟨ A⟩(t) = ⟨ A⟩_0 + δ⟨
        A⟩(t) + δ^2⟨
        A⟩(t)+…,
    
    where ⟨ A⟩_0 is the unperturbed average of A, and δ^n⟨
    A⟩(t) is the correction to the average of A to n-th order in the perturbation f. 
    The linear response δ⟨
    A⟩(t) can be expressed using the standard Kubo formula, which we re-derive for completeness in Appendix <ref>. 
    Going to one higher order in perturbation theory, we can write the second-order response δ^2⟨
    A⟩(t) in terms of a nonlinear response function
    χ^(2)_AB, defined via the constitutive relation
    
        δ^2⟨ A⟩(t)
           =∫ dt'dt”χ^(2)_AB(t-t',t-t”)f(t')f(t”),
    
    As derived in Appendix <ref>, we can take the Fourier transform of Eq. (<ref>) to find
    
        δ^2⟨ A⟩(ω)=∫  dω_1dω_2 δ(ω-ω_1-ω_2)χ^(2)_AB(ω_1,ω_2)f(ω_1)f(ω_2),
    
    where 
    
        χ^(2)_AB(ω_1,ω_2)=-∫_0^∞du∫_0^∞dv ⟨A(u+v)B(v)B(0)⟩_0e^iω_12^+ ue^iω_2^+ v,
    
    and where we have introduced the shorthand
    
        ω_1^+    = ω_1+iϵ, 
        ω_2^+    = ω_2+iϵ, 
        ω_12^+   =ω_1^++ω_2^+.
    
    
    Eq. (<ref>) is the second-order generalization of the Kubo formula. 
    Equivalent expressions appeared in the literature as early as Ref. <cit.>, and equivalent formulas have been widely applied more recently to the study of electromagnetic response functions <cit.>. 
    Here, we rewrite the nonlinear Kubo formula more generally, and in terms of time variables u and v that makes causality manifest: note that in Eq. (<ref>)A(u+v) is evaluated at a later time than B(v), which itself is evaluated at a later time than B(0). 
    
    Note, importantly, that the function χ^(2)_AB(ω_1,ω_2) is not itself a direct experimental observable. 
    This is because the perturbing fields on the right hand side of Eq. (<ref>) are symmetric under the interchange of integration variables ω_1↔ω_2. 
    The response function computed from an experimental measurement of δ^2⟨ A⟩(ω) is thus not χ^(2)_AB(ω_1,ω_2) but instead
    
        χ̅^(2)_AB(ω_1,ω_2)≡1/2(χ^(2)_AB(ω_1,ω_2)+χ^(2)_AB(ω_2,ω_1)).
    
    With this in mind, we will use Eq. (<ref>) for deriving general properties of the response function, and will explicitly symmetrize using Eq. (<ref>) before presenting physically-meaningful results. 
    
    As we discuss in Appendix <ref>, causality implies that both χ^(2)_AB(ω_1,ω_2) and χ̅^(2)_AB(ω_1,ω_2) are analytic functions in the upper half ω_1 and ω_2 planes. 
    As such, we can express χ^(2)_AB(ω_1,ω_2) as an integral over a spectral density. 
    By using the Fourier representation of the step function  [see Eq. (<ref>)], we can write
    
        χ^(2)_AB(ω_1,ω_2) = 1/π^2∫  dα dβχ^(2)”_AB(α,β)/(α-ω_12^+)(β-ω_2^+).
    
    The spectral density χ^(2)”_AB(α,β) is given as the Fourier transformed correlation function
    
        χ^(2)”_AB(α,β) = 1/4∫  du dv ⟨A(u+v)B(v)B(0)⟩_0e^i(α u + β v).
    
    Eq. (<ref>) is the first main result of this work. 
    Note that χ^(2)”_AB(α,β) is the Fourier transform of the correlation function appearing in the time-domain expression Eq. (<ref>) of the nonlinear response functions. 
    Thus, the spectral density also has the interpretation of an instantaneous response function, in analogy with the spectral density for linear response functions <cit.>. 
    Additionally, similar to the spectral density for linear response functions [see Eq. (<ref>) in Appendix <ref>], we can write χ^(2)”_AB(α,β) in the Lehmann representation as a weighted sum of Dirac delta functions. 
    We consider a ground state stationary density matrix of the form
    
        ρ_0 = ∑_n p_nnn,
    
    where |n⟩ is a complete set of eigenstates of the unperturbed Hamiltonian H_0 with energy E_n, each occupied with probability p_n. 
    Inserting complete sets of states in Eq. (<ref>) and evaluating the average we find
    
        χ^(2)”_AB(α,β) = π^2 ∑_nmℓ   [ℓAnnBmmBℓδ(α+E_ℓ-E_n)δ(β+E_ℓ-E_m)(p_ℓ-p_m). 
           +.ℓBnnBmmAℓδ(α+E_m-E_ℓ)δ(β+E_n-E_ℓ)(p_ℓ-p_n)].
    
    
    The spectral representation Eq. (<ref>) further allows us to derive generalized Kramers-Kronig relations for the second-order response function. 
    Using the Plemelj formula [see Eq. (<ref>)] to rewrite the denominators in Eq. (<ref>) and exploiting the analyticity of χ^(2)_AB(ω_1,ω_2) in the upper half ω_1 and ω_2 planes yields the decomposition
    
        χ^(2)_AB(ω_1,ω_2)    =  χ_AB^(2),1(ω_12,ω_2) + iχ_AB^(2),2(ω_12,ω_2),
    
    where we have defined 
    
        χ_AB^(2),1(x,y)   =1/π^2P∫  dα dβχ^(2)”_AB(α,β)/(α-x)(β-y)
        -χ^(2)”_AB(x,y),
        χ_AB^(2),2(x,y)   =1/πP∫  dα[χ^(2)”_AB(α,y)/α-x + χ^(2)”_AB(x,α)/α-y],
    
    and used P to denote the Cauchy principal value. 
    Symmetrizing Eqs. (<ref>) and (<ref>) using Eq. (<ref>), we can write for the physical response function
    
        χ̅^(2)_AB(ω_1,ω_2) = χ̅_AB^(2),1(ω_1,ω_2) + iχ̅_AB^(2),2(ω_1,ω_2)
    
    with
    
        χ̅_AB^(2),1(ω_1,ω_2)   =1/2( χ_AB^(2),1(ω_12,ω_1)+ χ_AB^(2),1(ω_12,ω_2) ), 
        χ̅_AB^(2),2(ω_1,ω_2)   =1/2( χ_AB^(2),2(ω_12,ω_1)+ χ_AB^(2),2(ω_12,ω_2) ).
    
    
    Unfortunately, unlike in the linear Kramers-Kronig relation, neither χ_AB^(2),1(ω_12,ω_2) nor χ_AB^(2),2(ω_12,ω_2) are proportional to the spectral density χ^(2)”_AB(ω_12,ω_2), meaning that direct measurement of the nonlinear spectral density χ^(2)”_AB(α,β) is not straightforward. 
    Nevertheless, χ_AB^(2),1(ω_12,ω_2) and χ_AB^(2),2(ω_12,ω_2) still form a Hilbert transform pair, generalizing the linear Kramers-Kronig relations. 
    In particular, 
    
        χ_AB^(2),1(x,y)   =1/πP∫ dαχ_AB^(2),2(α,y)/α-x
           = 1/πP∫ dαχ_AB^(2),2(x,α)/α-y, 
        χ_AB^(2),2(x,y)   =-1/πP∫ dαχ_AB^(2),1(α,y)/α-x
           = -1/πP∫ dαχ_AB^(2),1(x,α)/α-y.
    
    Similar generalized Kramers-Kronig relations have been previously obtained in, e.g. 
    Refs. <cit.>, although the connection with the spectral density was not emphasized in those works. 
    Note that, because these relationships do not directly involve the absorptive and reactive parts of the response, they cannot be applied to nonlinear spectroscopic data in the same, straightforward way that the usual Kramers-Kronig relations can be applied to linear spectroscopies. 
    However, they could be implemented in an iterative, self-consistent scheme to accomplish the same purpose.
    
    Although the second-order spectral density Eq. (<ref>) is not directly measurable through absorption, it can still be probed experimentally through the high-frequency behavior of the response function. 
    In particular, we will see now that the high-frequency asymptotic decay of χ̅^(2)(ω_1,ω_2) is determined by frequency moments of the spectral density χ^(2)”(α,β).
    
    
    
    § ASYMPTOTICS AND SUM RULES FOR SECOND HARMONIC GENERATION
    
    
    One main utility of the spectral representation Eq. (<ref>) is that it allows us to derive sum rules for the large frequency asymptotic behavior of χ̅^(2)(ω_1,ω_2). 
    Formally, we can expand Eq. (<ref>) in the limit of large ω_1 and ω_2 of the same sign and, after symmetrizing under the exchange of ω_1↔ω_2 [as per Eq. (<ref>)] we obtain the asymptotic expansion
    
        χ̅^(2)(ω_1→∞,ω_2→∞)∼1/2π^2∑_n,m=0^∞1/(ω_1+ω_2)^n+1(1/ω_1^m+1+1/ω_2^m+1)∫  dα dβ[α^nβ^m χ^(2)”_AB(α,β)].
    
    Since this asymptotic expansion only holds for both ω_1 and ω_2 large, its physical significance is somewhat obscure. 
    However, we can apply Eq. (<ref>) to the case of second harmonic generation ω_1=ω_2=ω. 
    We define the second harmonic generation response function
    
        χ_AB^SHG(ω) = χ̅^(2)(ω,ω) = χ^(2)(ω,ω)
    
    which determines the leading order contribution to the 2ω Fourier component of ⟨ A⟩ in response to a harmonic drive f(ω) at frequency ω via
    
        ⟨ A⟩(2ω) = χ_AB^SHG(ω)[f(ω)]^2 + 𝒪(f^4).
    
    Setting ω_1=ω_2=ω in the asymptotic expansion Eq. (<ref>), we find that the second harmonic generation response at large frequency is given by
    
        χ̅^SHG_AB(ω→∞)
           ∼∑_N=0^∞ν^(N)_AB/2ω^N+2,
    
    where we have defined the weighted moments ν^(N)_AB as
    
        ν^(N)_AB=1/π^2∑_m=0^N ∫ dα dβ(α/2)^N-mβ^mχ^(2)”_AB(α,β).
    
    We thus see that the moments ν^(N)_AB of the nonlinear spectral density χ^(2)”_AB(αβ) determine the asymptotic decay of the second harmonic generation response at large frequency. 
    Furthermore, since χ^(2)”_AB(α,β) is given in Eq. (<ref>) as the Fourier transform of a correlation function, we can express the moments ν^(N)_AB in terms of equal time correlation functions. 
    For concreteness, we can evaluate the first three moments N=0,1 and 2 explicitly. 
    From Eqs. (<ref>) and (<ref>) we find
    
        ν^(0)_AB   = 1/π^2∫ dα dβχ^(2)”_AB(αβ) 
        
           =⟨ABB⟩_0.
    
    Similarly, for ν^(1)_AB we find
    
        ν^(1)_AB   = 1/π^2∫ dα dβ(α/2+β)χ^(2)”_AB(α,β)
        
           =⟨i∂_tABB⟩_0 + 1/2⟨i∂_t ABB⟩_0.
    
    Finally, for ν^(2)_AB we have
    
        ν^(2)_AB   = 1/π^2∫ dα dβ(α^2/4+β^2+αβ/2)χ^(2)”_AB(αβ) 
           =-⟨∂^2_tABB⟩_0 - 1/2⟨∂_t∂_t ABB⟩_0 
           -1/4⟨∂^2_t ABB⟩_0.
    
    
    Eqs. (<ref>)–(<ref>) relate the frequency moments of the spectral density χ^(2)”_AB to ground state averages of commutators of A, B, and their time derivatives. 
    They are nonlinear generalizations of the standard sum rules for the linear response coefficients introduced in Refs. <cit.>, and reviewed in Appendix <ref>. 
    Comparing with the time-domain expression for the nonlinear response function [i.e., Eq. (<ref>) and the integrand in Eq. (<ref>)], we see that the averages appearing in the sum rule correspond to a gradient expansion of the time-domain response near u=v=0, sometimes called the instantaneous response function <cit.>. 
    
    We see that the sum rules constrain the large frequency behavior of the second harmonic generation response via Eq. (<ref>), which can be measured via nonlinear spectroscopic techniques. 
    Furthermore, the sum rules and asymptotic expansions derived here make no assumptions on the unperturbed Hamiltonian other than that it conserves energy, and that averages are taken in a state described by a stationary density matrix. 
    In particular, the asymptotic expansion Eq. (<ref>) and the sum rules Eqs. (<ref>)–(<ref>) (as well as their generalization to higher moments) hold for arbitrary interacting many-body systems. 
    In the next section, we will apply this formalism to analyze the response of the electronic charge density ρ_𝐪 to the scalar electric potential ϕ_𝐪(t).
    
    
    
    § DENSITY RESPONSE AND THE SECOND-ORDER F-SUM RULE
    
    
    We will now study the implications of our nonlinear sum rules Eqs. (<ref>)–(<ref>) for second-order density response. 
    We consider the response of the charge density operator
    
        A=ρ_𝐪
    
    to a scalar electric potential perturbation
    
        H_1=∑_𝐪ϕ_𝐪(t)ρ_-𝐪,
    
    where ϕ_𝐪(t) is that (spatial) Fourier transform of the electric potential. 
    Following our derivation of the nonlinear response function in Appendix <ref>, we can define the second order density-density response function χ^(2)_ρρ(ω_1,𝐪_1,ω_2,𝐪_2) via
    
        δ^2⟨ρ_𝐪⟩(ω) = ∑_𝐪_1,𝐪_2∫ dω_1 dω_2χ^(2)_ρρ(ω_1,𝐪_1,ω_2,𝐪_2)δ(𝐪-𝐪_12)δ(ω-ω_12)f_𝐪_1(ω_1)f_𝐪_2(ω_2),
    
    where the nonlinear Kubo formula gives
    
        χ^(2)_ρρ(ω_1,𝐪_1,ω_2,𝐪_2)= -∫_0^∞du∫_0^∞dv ⟨ρ_𝐪_12(u+v)ρ_-𝐪_1(v)ρ_-𝐪_2(0)⟩_0e^iω_12^+ ue^iω_2^+ v,
    
    and we have introduced 
    
        𝐪_12=𝐪_1+𝐪_2
    
    by analogy with Eq. (<ref>). 
    As in the general case Eq. (<ref>), the experimentally-measurable response function is not directly χ^(2)_ρρ(ω_1,𝐪_1,ω_2,𝐪_2), but is instead the symmetrized
    
        χ̅^(2)_ρρ(ω_1,𝐪_1,ω_2,𝐪_2)    = 1/2(χ^(2)_ρρ(ω_1,𝐪_1,ω_2,𝐪_2)+χ^(2)_ρρ(ω_2,𝐪_2,ω_1,𝐪_1)).
    
    Following the logic of Eq. (<ref>), we can introduce a spectral representation
    
        χ^(2)_ρρ(ω_1,𝐪_1,ω_2,𝐪_2) = 1/π^2∫  dα dβχ^(2)”_ρρ(α,𝐪_1,β,𝐪_2)/(α-ω_12^+)(β-ω_2^+),
    
    with the spectral density χ^(2)”_ρρ(α,𝐪_1,β,𝐪_2) defined as the Fourier transform
    
        χ^(2)”_ρρ(α,𝐪_1,β,𝐪_2) = 1/4∫  du dv ⟨ρ_𝐪_12(u+v)ρ_-𝐪_1(v)ρ_-𝐪_2(0)⟩_0e^i(α u + β v).
    
    Let us now consider second harmonic generation response ω_1=ω_2=ω to apply our sum rules Eqs. (<ref>)–(<ref>). 
    To begin, we can define the second harmonic generation response function
    
        χ̅^SHG_ρρ(ω,𝐪_1,𝐪_2)    = χ̅^(2)_ρρ(ω,𝐪_1,ω,𝐪_2) 
           = 1/2(χ^(2)_ρρ(ω,𝐪_1,ω,𝐪_2).
           .+χ^(2)_ρρ(ω,𝐪_2,ω,𝐪_1)).
    
    We can expand the spectral density representation Eq. (<ref>) to obtain an asymptotic series for χ̅^SHG_ρρ(ω→∞,𝐪_1,𝐪_2) in the limit of large ω. 
    In analogy with Eq. (<ref>) for the general response function, we find
    
        χ̅^SHG_ρρ(ω→∞,𝐪_1,𝐪_2)∼∑_N=0^∞ν_ρρ,𝐪_1,𝐪_2^(N)+ν_ρρ,𝐪_2,𝐪_1^(N)/4ω^N+2,
    
    where the moments ν_ρρ,𝐪_1,𝐪_2^(N) are defined via
    
        ν_ρρ,𝐪_1,𝐪_2^(N) = 1/π^2∑_m=0^N∫ dα dβ   (α/2)^N-mβ^m 
           ×χ^(2)”_ρρ(α,𝐪_1,β,𝐪_2).
    
     For the zeroth moment ν^(0)_ρρ,𝐪_1,𝐪_2 we find
    
        ν^(0)_ρρ,𝐪_1,𝐪_2    =  1/π^2∫ dα dβχ^(2)”_ρ,ρ(α,𝐪_1,β,𝐪_2) 
           =⟨ρ_𝐪_12ρ_-𝐪_1ρ_-𝐪_2⟩_0 
           =0,
    
    since all Fourier components of the density operator are mutually commuting. 
    For the first moment ν^(1)_ρρ, 𝐪_1,𝐪_2 we find from Eq. (<ref>) that
    
        ν^(1)_ρρ, 𝐪_1,𝐪_2    =  1/π^2∫ dα dβ(α/2+β)χ^(2)”_ρ,ρ(α,𝐪_1,β,𝐪_2) 
           =1/2⟨i∂_t ρ_𝐪_12ρ_-𝐪_1ρ_-𝐪_2⟩_0,
    
    where we have again made use of the fact that density operators have vanishing equal-time commutators. 
    
    We can additionally use charge conservation to re-express Eq. (<ref>) in terms of a variation of the unperturbed Hamiltonian H_0 with respect to the electromagnetic vector potential. 
    As we review in Appendix <ref>, conservation of charge along with the canonical commutation relations imply that for any operator 𝒪,
    
        [𝒪,ρ_𝐪] =.q^μδ𝒪_A/δ A^μ_-𝐪|_𝐀→ 0,
    
    where 𝒪_A is the operator 𝒪 minimally coupled to the external vector potential 𝐀_𝐪.
    We find
    
        ν^(1)_ρρ, 𝐪_1,𝐪_2   = 1/2(q_12^μ)(q_1^ν)(q_2^λ)
           ×.⟨δ^3 H/δ A^μ_-𝐪_12δ A^ν_𝐪_1δ A^λ_𝐪_2⟩_0|_𝐀→ 0,
    
    The third variation of the Hamiltonian with respect to the vector potential gives the generalized diamagnetic current, describing corrections to the current operator that are quadratic in the vector potential. 
    These are sometimes known as “three photon vertices” <cit.>. 
    As shown in Ref. <cit.>, the right hand side of Eq. (<ref>) can be evaluated as a nested commutator between the position operator and the ordinary (paramagnetic) current operator.
    
    
    Combining Eqs. (<ref>) and (<ref>), we see that the second harmonic generation response function satisfies
    
        lim_ω→∞ω^3   χ̅^SHG_ρρ(ω,𝐪_1,𝐪_2) =  
           1/4(q_12^μ)(q_1^ν)(q_2^λ).⟨δ^3 H/δ A^μ_-𝐪_12δ A^ν_𝐪_1δ A^λ_𝐪_2⟩_0|_𝐀→ 0.
    
    Eqs. (<ref>) and (<ref>) give the first generalization of the linear f-sum rule [see Eq. (<ref>)] to the second-order density response, and are a main result of this note. 
    
    
    However, for nonrelativistic systems, the Hamiltonian is at most quadratic in momentum and interactions are momentum-independent. 
    We can write a general nonrelativistic Hamiltonian as
    
        H_n.r.=∑_i |𝐩_i|^2/2m - λ⃗(⃗𝐱⃗_⃗i⃗)⃗·(𝐩×σ⃗) + V(𝐱_i) + 1/2∑_i≠ j U(𝐱_i-𝐱_j),
    
    where i,j index particles in the system, m is the electron mass, V is an external potential, λ⃗ is the spin-orbit potential, σ⃗ is a vector of Pauli matrices acting on the electron spins, and U is a pair interaction. 
    Minimally coupling H_𝐧.𝐫. to a vector potential via 𝐩_i→𝐩_i - 𝐀(𝐱_i) produces terms of at-most quadratic order in the vector potential 𝐀. 
    As such, we have
    
        .δ^3 H_n.r./δ A^μ_-𝐪_12δ A^ν_𝐪_1δ A^λ_𝐪_2|_𝐀→ 0 = 0.
    
    Thus, for systems governed by nonrelativistic Hamiltonians of the form of Eq. (<ref>), ν^(1)_ρρ, 𝐪_1,𝐪_2 vanishes identically, and we need to go to the second moment to obtain a nontrivial sum rule and nonzero asymptotic expansion. 
    
    Turning to ν^(2)_ρρ,𝐪_1,𝐪_2 we have
    
        ν^(2)_ρρ,𝐪_1,𝐪_2 = 1/π^2∫ dα dβ(α^2/4+β^2+αβ/2)χ^(2)”_ρρ(α,𝐪_1,β,𝐪_2) 
           = -2⟨∂_t∂_t ρ_𝐪_12ρ_-𝐪_1ρ_-𝐪_2⟩_0 +⟨∂^2_t ρ_𝐪_12ρ_-𝐪_1ρ_-𝐪_2⟩_0/4
           = -3⟨∂_t∂_t ρ_𝐪_12ρ_-𝐪_1ρ_-𝐪_2⟩_0 +⟨∂_t ρ_𝐪_12∂_tρ_-𝐪_1ρ_-𝐪_2⟩_0/4.
    
    Eq. (<ref>) is general, and requires no assumptions on the form of the Hamiltonian. 
    However, for nonrelativistic systems we can simplify it further by making use of Eqs. (<ref>) and Eq. (<ref>) to find that for nonrelativistic systems the density algebra satisfies
    
        ∂_t ρ_𝐪ρ_𝐪'→_n.r. i𝐪·𝐪'/mρ_𝐪+𝐪'.
    
    By using the Jacobi identity to rewrite the second term in Eq. (<ref>) solely in terms of commutators of the density with time derivatives of the density, we find that in the nonrelativistic limit
    
        ν^(2)_ρρ,𝐪_1,𝐪_2   →_n.r.3 i𝐪_12·𝐪_1 /4m⟨∂_t ρ_𝐪_2ρ_-𝐪_2⟩_0
            -1/4⟨∂_tρ_-𝐪_1ρ_-𝐪_2∂_tρ_𝐪_12⟩_0 
           +1/4⟨∂_tρ_𝐪_12ρ_-𝐪_2∂_tρ_-𝐪_1⟩_0 
           =n̅𝐪_12·[3 𝐪_1 |𝐪_2|^2 +𝐪_12(𝐪_1·𝐪_2) + 𝐪_2 |𝐪_1|^2] /4m^2,
    
    where n̅=⟨ρ_0⟩_0 is the average charge density in equilibrium.
    
    Substituting Eq. (<ref>) into our asymptotic expansion Eq. (<ref>), we see that, in the nonrelativistic limit the second harmonic generation response at large frequency is asymptotically given by
    
        χ̅^SHG_ρρ(ω→∞,𝐪_1,𝐪_2)   ∼_n.r.n̅/8m^2ω^4(4|𝐪_1|^2|𝐪_2|^2.
           .-4(𝐪_1·𝐪_2)^2.
           .+3𝐪_1·𝐪_2|𝐪_12|^2) +𝒪(ω^-5).
    
    Finally, note that for a harmonic perturbation, f_𝐪 has only two nonvanishing Fourier components at 𝐪=±𝐪_0. 
    In this case, the relevant moments are
    
        ν^(2)_ρρ(𝐪_0,-𝐪_0) →_n.r. 0
    
    and
    
        ν^(2)_ρρ(𝐪_0,𝐪_0) →_n.r.3|𝐪_0|^4n̅/m^2
    
    and hence
    
        χ̅^SHG_ρρ(ω→∞,𝐪_0,-𝐪_0)   ∼_n.r. 0, 
        χ̅^SHG_ρρ(ω→∞,𝐪_0,𝐪_0)   ∼_n.r.3n̅|𝐪_0|^4/2m^2ω^4,
    
    Note that χ̅^SHG_ρρ(ω,𝐪_0,-𝐪_0) determines the second order response 
    
        δ^2⟨ρ_0⟩(2ω) = χ̅^SHG_ρρ(ω,𝐪_0,-𝐪_0)f_𝐪_0(ω)f_-𝐪_0(ω).
    
    In fact, χ̅^SHG_ρρ(ω,𝐪_0,-𝐪_0) is identically zero for systems with particle number conservation, since in that case the 𝐪=0 Fourier component of the density operator is time independent, so the commutator in Eq. (<ref>) vanishes.
    
    Similarly, χ̅^SHG_ρρ(ω,𝐪_0,𝐪_0) determines the second-order response
    
        δ^2⟨ρ_2𝐪_0⟩(2ω) = χ̅^SHG_ρρ(ω,𝐪_0,𝐪_0)f_𝐪_0(ω)f_𝐪_0(ω).
    
    Thus, the nonrelativistic sum rules Eqs. (<ref>) and (<ref>) constrain the large-frequency asymptotic response of the density Fourier component ⟨ρ_2𝐪_0⟩(2ω) in the limit of large frequency. 
    These sum rules place important constraints on, and should aid the interpretation of, the newest generation of nonlinear x-ray optics experiments, which exploit second-order optical effects at large momenta and have potential new applications in surface science, low-dimensional transport, and multidimensional spectroscopies involving core levels <cit.>. 
    
    
    
    § EXAMPLE: DENSITY RESPONSE FOR FREE FERMIONS IN A CRYSTAL
    
    
    We will now see how our density response sum rules constrain the properties of noninteracting electrons in solids. 
    We consider an unperturbed Hamiltonian H_0 with discrete translation symmetry consisting solely of one-particle terms. 
    We can label the single particle eigenstates |ψ_n𝐤⟩ of H_0 by their crystal momentum 𝐤 in the first Brillouin zone, such that
    
        H_0|ψ_n𝐤⟩ = ϵ_n𝐤|ψ_n𝐤⟩.
    
    It will be computationally convenient to recast Eq. (<ref>) in terms of the cell-periodic functions
    
        |u_n𝐤⟩ = e^-i𝐤·𝐱|ψ_n𝐤⟩,
    
    where 𝐱 is the single-particle position operator. 
    Introducing the Bloch Hamiltonian
    
        H_𝐤 = e^-i𝐤·𝐱H_0e^i𝐤·𝐱,
    
    we have
    
        H_𝐤|u_n𝐤⟩ = ϵ_n𝐤|u_n𝐤⟩.
    
    We will normalize the |u_n𝐤⟩ over a single unit cell. 
    Introducing creation and annihilation operators c^†_n𝐤, c_n𝐤 corresponding to the states |ψ_n𝐤⟩, we can write the Hamiltonian H_0 in second quantization as
    
        H_0 = ∑_n𝐤ϵ_n𝐤c^†_n𝐤c_n𝐤.
    
    Similarly, the density operator in second quantization is given by
    
        ρ_𝐪 = ∑_nm𝐤⟨u_n𝐤|u_m𝐤+𝐪⟩c^†_n𝐤c_m𝐤+𝐪,
    
    where the overlap between cell-periodic functions is defined as an integral over one unit cell. 
    
    For a free-fermion ground state, we can directly evaluate the spectral density Eq. (<ref>) for the second-order density response function as a sum of density operator matrix elements times Fermi-Dirac distribution functions, using Eq. (<ref>). 
    The expression is exact, but unwieldy. 
    On the other hand, we can use Eqs. (<ref>) and (<ref>) along with canonical anticommutation relations to derive a simple exact expression for the nonlinear sum rule Eq. (<ref>). 
    To do so, we first use the equations of motion to rewrite Eq. (<ref>) as
    
        ν^(1)_ρρ,𝐪_1,𝐪_2 = -1/2Hρ_𝐪_12ρ_-𝐪_1ρ_-𝐪_2_0.
    
    To evaluate the nested commutators, we can make use of the following relation: Given any single-body operator
    
        B_𝐪=∑_𝐤nmu_n𝐤B_𝐤,𝐪u_m𝐤+𝐪c^†_n𝐤c_n𝐤+𝐪,
    
    we can use the canonical anticommutation relations to find
    
        [B_𝐪,ρ_𝐪'] = ∑_𝐤nmu_n𝐤B_𝐤,𝐪-B_𝐤+𝐪',𝐪u_m𝐤+𝐪+𝐪'.
    
    Making extensive use of Eq. (<ref>) to evaluate the nested commutators in Eq. (<ref>), we find that
    
        Hρ_𝐪_12ρ_-𝐪_1ρ_-𝐪_2    = ∑_nm𝐤c^†_n𝐤c_m𝐤**u_n𝐤H_𝐤-𝐪_1-𝐪_2 + H_𝐤+𝐪_1+H_𝐤+𝐪_2-H_𝐤+𝐪_1+𝐪_2-H_𝐤-𝐪_1-H_𝐤-𝐪_1u_m𝐤.
    
    
    We can evaluate this average in a thermal ground state with fixed temperature T and chemical potential μ. 
    Let f_n𝐤 be the Fermi-Dirac distribution function evaluated at energy ϵ_n𝐤. 
    Then we have
    
        ν^(1)_ρρ,𝐪_1,𝐪_2   = 1/2∑_n𝐤f_n𝐤**u_n𝐤H_𝐤+𝐪_1+𝐪_2 + H_𝐤-𝐪_1+H_𝐤-𝐪_2-H_𝐤-𝐪_1-𝐪_2-H_𝐤+𝐪_1-H_𝐤+𝐪_1u_n𝐤.
    
    Eq. (<ref>) expresses the first moment of the second order spectral density entirely in terms of matrix elements of the Bloch Hamiltonian. 
    For low-energy Wannier-based models of condensed matter systems, Eq. (<ref>) combined with Eq. (<ref>) shows that there is an effective second-order diamagnetic current proportional to ν^(1)_ρρ,𝐪_1,𝐪_2. 
    Furthermore, we can use the Karplus-Schwinger relation of Refs. <cit.> to rewrite Eq. (<ref>) in terms of derivatives of the Bloch Hamiltonian as
    
        ν^(1)_ρρ,𝐪_1,𝐪_2 =    q_12^μ q_1^ν q_2^λ1/2[∑_n𝐤∫_0^1dλ∫_0^1dλ_1∫_0^1dλ_2 f_n𝐤.
           .×**u_n𝐤∂^3 H_𝐤-λ𝐪_12+λ_1𝐪_1+λ_2𝐪_2/∂ k^μ∂ k^ν∂ k^λu_n𝐤],
    
    which shows that the longitudinal component of the diamagnetic current coincides with the longitudinal component of the term in brackets in Eq. (<ref>). 
    This is consistent with the results of Ref. <cit.>.
    
    
    
    § SECOND-ORDER CONDUCTIVITY SUM RULE
    
    
    From charge conservation, it follows that the second order density response function χ^(2)_ρρ(𝐪_1,ω_1,𝐪_2,ω_2) determines the longitudinal components of the second order conductivity tensor, and hence our sum rules in Eqs. (<ref>) and (<ref>)–(<ref>) for the density response imply longitudinal second order conductivity sum rules. 
    To see this concretely, let us return to the general formalism of Sec. <ref>. 
    We consider again a scalar potential perturbation
    
        H_1=∑_𝐪ϕ_𝐪(t)ρ_-𝐪,
    
    and look at the response of the current operator 𝐣_𝐪. 
    At second order in the scalar potential, we have
    
        δ^2⟨ j^μ_𝐪⟩(ω) = ∑_𝐪_1,𝐪_2∫ dω_1 dω_2χ^(2),μ_jρ(ω_1,𝐪_1,ω_2,𝐪_2)δ(𝐪-𝐪_12)δ(ω-ω_12)ϕ_𝐪_1(ω_1)ϕ_𝐪_2(ω_2),
    
    where the nonlinear Kubo formula gives
    
        χ^(2),μ_jρ(ω_1,𝐪_1,ω_2,𝐪_2)= -∫_0^∞du∫_0^∞dv ⟨j^μ_𝐪_12(u+v)ρ_-𝐪_1(v)ρ_-𝐪_2(0)⟩_0e^iω_12^+ ue^iω_2^+ v,
    
    However, from gauge invariance, we know that the current cannot respond to the scalar potential directly, but only to the electric field. 
    In the gauge of Eq. (<ref>) where the vector potential is zero, the electric field is given by
    
        𝐄_𝐪(t) = -i𝐪ϕ_𝐪(t).
    
    Defining the nonlinear conductivity σ^μνλ(ω_1,𝐪_1,ω_2,𝐪_2) via the constitutive relation
    
        δ^2⟨ j^μ_𝐪⟩(ω)=∑_𝐪,𝐪'∫ dω_1 dω_2σ^μνλ(ω_1,𝐪_1,ω_2,𝐪_2)δ(𝐪-𝐪_12)δ(ω-ω_12)E^ν_𝐪_1(ω_1)E^λ_𝐪_2(ω_2),
    
    We see that we can choose [There is an ambiguity here due to the fact that only the response functions symmetrized under 𝐪_1↔𝐪_2 are determined from the constitutive relations. Nevertheless, we are free to make any choice we wish for the antisymmetric terms provided we symmetrize before obtaining physical results.]
        χ^(2),μ_jρ(ω_1,𝐪_1,ω_2,𝐪_2)=-q_1^ν q_2^λσ^μνλ(ω_1,𝐪_1,ω_2,𝐪_2).
    
    To progress further, we can contract χ^(2),μ_jρ with 𝐪_12 and use the continuity equation
    
        i∂_tρ_𝐪_12 = 𝐪_12·𝐣_𝐪_12
    
    to relate σ^μνλ(ω_1,𝐪_1,ω_2,𝐪_2) with the density response function χ^(2)_ρρ(ω_1,𝐪_1,ω_2,𝐪_2) from Eq. (<ref>). 
    Inserting Eq. (<ref>) into Eq. (<ref>) and integrating by parts we find
    
        -q_12^μ q_1^ν q_2^λσ^μνλ(ω_1,𝐪_1,ω_2,𝐪_2) = ω_12^+χ^(2)_ρρ(ω_1,𝐪_1,ω_2,𝐪_2),
    
    Eq. (<ref>) relates the longitudinal second-order conductivity to the second-order density response function. 
    Going further, we can make use of the spectral representation Eq. (<ref>) to express the longitudinal second-order conductivity as 
    
        -q_12^μ q_1^ν q_2^λ   σ^μνλ(ω_1,𝐪_1,ω_2,𝐪_2) 
           = 1/π^2∫  dα dβαχ^(2)”_ρρ(α,𝐪_1,β,𝐪_2)/(α-ω_12^+)(β-ω_2^+)
           =-1/π^2∫  dα dβσ_L”(α,𝐪_1,β,𝐪_2)/(α-ω_12^+)(β-ω_2^+).
    
    We see then that the spectral density σ_L”(α,𝐪_1,β,𝐪_2) for the longitudinal second-order conductivity is given by -α times the spectral density χ^(2)”_ρρ(α,𝐪_1,β,𝐪_2) for the second-order density response function. 
    We can now follow the logic of Sec. <ref> to asymptotically expand Eq. (<ref>) at large frequency to deduce sum rules for the longitudinal conductivity. 
    Focusing on the second harmonic current ω_1=ω_2=ω at large ω and explicitly symmetrizing under 𝐪_1↔𝐪_2 to obtain the physical response [we denote the symmetrization with an overbar, extending the notation of Eq. (<ref>)]
    
        q_12^μ q_1^ν q_2^λσ̅^μνλ(ω,𝐪_1,ω,𝐪_2) ∼∑_N=0^∞μ^(N)_𝐪_1,𝐪_2+μ^(N)_𝐪_2,𝐪_1/4ω^N+2,
    
    where the moments μ^(N)_𝐪_1,𝐪_2 are defined as
    
        μ^(N)_𝐪_1,𝐪_2   =1/π^2∑_m=0^N∫ dα dβ(α/2)^N-mβ^mσ_L”(α,𝐪_1,β,𝐪_2) 
           =1/π^2∑_m=0^N∫ dα dβ(α/2)^N-mβ^m
        αχ^(2)”_ρρ(α,𝐪_1,β,𝐪_2).
    
    In particular, we find that
    
        μ^(0)_𝐪_1,𝐪_2   =2ν^(1)_ρρ, 𝐪_1,𝐪_2, 
        μ^(1)_𝐪_1,𝐪_2   =2ν^(2)_ρρ,𝐪_1,𝐪_2,
    
    where ν^(1)_ρρ, 𝐪_1,𝐪_2 and ν^(2)_ρρ,𝐪_1,𝐪_2 are the moments defined in Eqs. (<ref>) and (<ref>), respectively. 
    Eqs. (<ref>) and (<ref>) are sum rules for the longitudinal second order conductivity. 
    In particular, Eq. (<ref>) shows that μ^(0)_𝐪_1𝐪_2 is proportional to the average third derivative of the Hamiltonian with respect to the vector potential. 
    This is equivalent to the instantaneous sum rule for the second-order conductivity derived in Ref. <cit.>. 
    Our formalism extends this result beyond leading asymptotic order in ω, which is particularly important for nonrelativistic systems in which ν^(1)_ρρ, 𝐪_1,𝐪_2 is identically zero. 
    Finally, applied to free-fermion systems, we see from Eqs. (<ref>) and (<ref>) that the conductivity sum rule Eq. (<ref>) can be written entirely in terms of matrix elements of the Bloch Hamiltonian (or of its third derivative).
    
    
    
    § RECTIFICATION SUM RULE
    
    Up to now, we have focused primarily on second-harmonic generation response ω_1=ω_2. 
    However, for response to real harmonic perturbations at second order, there is also a rectification response with ω_1=-ω_2. 
    The rectification response determines the zero frequency (time-independent) response of the average ⟨ A ⟩(ω=0) to the perturbation f(t) at second order; setting ω=0 in Eq. (<ref>), we can defined the rectification response χ_AB^R(ω) via
    
        ⟨ A ⟩(ω=0) = ∫ dωχ_AB^R(ω)|f(ω)|^2,
    
    where we have used the fact that for real perturbations f(-ω)=f(ω)^*. 
    Comparing with Eq. (<ref>) we have that
    
        χ_AB^R(ω) = 1/2(
        χ_AB^(2)(ω,-ω) + χ_AB^(2)(-ω,ω)
        ).
    
    The spectral representation for the rectification response takes a particularly simple form. 
    Setting ω_1=-ω_2=ω in the spectral representation Eq. (<ref>) and symmetrizing, we find
    
        χ_AB^R(ω)= 1/π^2∫ dα dββχ_AB^(2)”(α,β)/(α-2iϵ)(β^2-(ω^+)^2),
    
    which defines the rectification response in terms of an integral of the spectral density.
    
    Analogous to our treatment second-harmonic generation response in Sec. <ref>, we can use Eq. (<ref>) to derive an asymptotic expansion for the rectification response in the limit of large frequency |ω|→∞. 
    We find asymptotically that
    
        χ_AB^R(ω→∞)∼∑_N=0^∞ζ^(N)_AB/2ω^2N+2,
    
    where we have introduced the moments
    
        ζ^(N)_AB = -1/π^2∫ dα dββ^2N+1χ_AB^(2)”(α,β)/α-2iϵ.
    
    Note that, unlike our sum rule for second harmonic generation, the moments ζ^(N)_AB involve the negative first moment in α of the spectral density. 
    This is reminiscent of the sum rule for polarization fluctuations introduced in Refs. <cit.>. 
    In general, the presence of a negative power of α hinders our ability to simplify Eq. (<ref>) using our definition Eq. (<ref>) of the spectral density. 
    However, in the special case that A can be written as a time derivative
    
        A=i[H_0,C]≡∂_t C
    
    of some operator C in the Heisenberg picture, we can combine Eqs. (<ref>), (<ref>) and Eq. (<ref>) and integrate by parts to obtain
    
        ζ^N_AB = (-1)^NCB∂_t^2N+1B_0,
    
    where the time derivatives of B are evaluated in the Heisenberg picture as commutators with H_0. 
    Eqs. (<ref>), (<ref>), and (<ref>) give a family of generalized sum rules for the second order rectification response. 
    While not directly applicable to density response (since ρ_𝐪 cannot in general be written as the time derivative of an operator), our sum rule can be applied to the spatially uniform (𝐪→ 0) conductivity, using the fact that the uniform current can be written as the time derivative of the position operator. 
    Note that, unlike the rectification sum rule derived in Ref. <cit.>, the integral over α in Eq. (<ref>) extends over the entire real line rather than just the positive interval [0,∞). 
    While we can use our explicit form of the spectral density in Eq. (<ref>) to compute moments of the spectral density integrated over only positive frequencies, these are not directly related to the high-frequency asymptotics of the response function. 
    We will defer treatment of sum rules of that type to future work.
    
    
    
    
    § OUTLOOK
    
    
    In this work, we have examined the general theory of second-order response functions. 
    First, in Sec. <ref>, we exploited causality to introduce a spectral density representation of the second-order Kubo formula, generalizing the familiar Kramers-Kronig relations from linear response theory. 
    The spectral representation Eq. (<ref>) allowed us to express the large-frequency asymptotics of the second-order response function in terms of frequency moments of the spectral density. 
    To relate these moments to experimentally-observable responses, we focused on second-harmonic generation. 
    By using our explicit expression Eq. (<ref>) for the spectral density to evaluate the frequency moments, we showed how they could be expressed in terms of averages of equal-time commutators. 
    These relations are generalized sum rules relating equal-time averages, frequency moments of the spectral density, and large-frequency asymptotics of the second-harmonic generation response. 
    While we give explicit expressions for the first three sum rules in Eqs. (<ref>)–(<ref>), our method in Sec. <ref> straightforwardly extends to an infinite family of sum rules. 
    Next, we specialize in Sec. <ref> to the experimentally-relevant case of second-order density response. 
    Applying our formalism to this case allowed us to derive nonlinear generalizations of the linear f-sum rule. 
    In particular, we found in Eq. (<ref>) that the leading-order behavior of the second-harmonic density response at large frequency is determined by the first moment of the spectral density; this can be expressed via our sum rule as the average second-order diamagnetic current (third derivative of the Hamiltonian with respect to vector potential). 
    For nonrelativistic systems, this is identically zero and we have to go to next-order to find that our nonlinear sum rule relates the second moment of the spectral density response to the average charge density in equilibrium, in analogy with the linear f-sum rule. 
    We applied our results to a system of noninteracting electrons in a periodic potential in Sec. (<ref>), showing that the leading-order sum rule can be expressed entirely in terms of matrix elements of the Bloch Hamiltonian. 
    Using charge conservation, we showed how our second-order density response sum rules imply an infinite family of sum rules for the longitudinal second order conductivity, generalizing the results of Ref. <cit.> to nonzero wavevector and to higher moments of the conductivity.
    Finally, we examined the rectification response, the zero frequency response to a harmonic perturbation at second order. 
    We showed that for observables that can be written as time derivatives of other operators, the rectification response at large frequencies is determined by a set of sum rules expressible as averages of equal-time commutators.  
    
    Our work has implications for future theoretical and experimental studies on quantum materials. 
    First, we re-emphasize that our spectral density and sum rules in Secs. <ref>, <ref>, <ref>, <ref>, and <ref> are valid for any (energy-conserving) quantum system, and require no assumptions on the strength of electron-electron interactions or on the form of the unperturbed Hamiltonian. 
    As such, our sum rules are formally exact, making them relevant to experiments on strongly-correlated as well as band materials. 
    The second-order density response at large frequencies is directly accessible via nonlinear X-ray scattering, placing important constraints on such measurements, and enabling experimenters to calibrate the absolute scale of nonlinear cross sections in the same manner as linear spectroscopies <cit.>. 
    Additionally, Refs. <cit.> showed that N-point density-density correlation functions in a Fermi liquid are related to multipartite entanglement and Fermi surface topology in the ground state. 
    Since these correlation functions appear in the Fourier transform of χ^(2)”_ρρ in the time-domain, the connection between these works and our sum rules is a fruitful area for future theoretical exploration.
    
    Furthermore our general method in Secs. <ref> and <ref> can be applied to derive sum rules for the conductivity tensor directly, rather than just its longitudinal components. 
    For optical response at zero wavevector, our formalism of Sec. <ref> can be applied directly using 𝐣_0 = i[H_0,𝐗], where 𝐗 is the many-body position operator, provided care is taken to also include the generalized diamagnetic conductivity (see, e.g. the generalized current vertices in Refs. <cit.>). 
    In this way our approach can bridge the gap between the Riemannian-geometric approach to nonlinear optical response from Ref. <cit.> and sum-rule based constraints generalizing the “quantum weight” of Refs. <cit.>.
    
    Finally, let us remark on the applicability of the sum rules presented here to low-energy effective models, where very high-energy degrees of freedom have been projected out. 
    In the context of such low-energy effective models, we typically consider the response of operators projected into the low-energy subspace (i.e. projected density operators), and restrict frequency integrations in sum rules to the range of some gap or cutoff scale. 
    For effective models whose low-energy subspace is spanned by ultralocalized Wannier functions, we expect our response functions and sum rules to hold in the form presented in Secs. <ref> and <ref>. 
    This is because the algebra of projected density operators is unmodified provided all moments of the position operator in the Wannier basis are diagonal, as is the case in the ultralocalized (or tight-binding) limit <cit.>. 
    However, in moiré and other strongly-interacting topological systems, we are often interested in projecting onto low-energy degrees of freedom with delocalized Wannier functions, which is a situation in which low-energy partial sum rules are known to break down <cit.>. 
    In such cases, a full theory of projected response functions is required. 
    Preliminary work in this direction has focused on the density algebra in fractional quantum Hall <cit.> and Chern insulators <cit.>. 
    Additionally, recent work on linear sum rules in Ref. <cit.> has shown the importance of the projected density operator on optical sum rules over restricted frequency range. 
    Using the tools we developed here, the extension of these results to nonlinear response functions can be undertaken.
    
    
    
    The authors thank C. L. 
    Kane, J. E. 
    Moore and P. 
    Phillips for fruitful discussions. 
    This work is supported by the U.S. 
    DOE, Office of Basic Energy Sciences, Energy Frontier Research Center for Quantum Sensing and Quantum Materials through Grant No. 
    DE-SC0021238. 
    B. B. received additional support for theoretical development from the Alfred P. 
    Sloan foundation, and the National Science Foundation under grant DMR-1945058. 
    P. A. gratefully acknowledges additional support from the EPiQS program of the Gordon and Betty Moore Foundation, grant no. 
    GBMF9452.
    
    
    § REVIEW OF LINEAR SUM RULES
    
    In this Appendix, we review the derivation of the spectral density and general sum rules for linear response functions. 
    We follow the approach of Refs. <cit.>. 
    Suppose our system of interest is governed by Hamiltonian H_0, and we apply a time-dependent external field that couples to an operator B via the perturbing Hamiltonian
    
        H_1 = f(t)B.
    
    We can then write the average of an operator A to linear order in B as
    
        δ⟨ A⟩(t) = ⟨ A ⟩ (t) -⟨ A⟩_0 = ∫ dt' χ_AB(t-t') f(t'),
     
    where ⟨ A⟩_0=tr(Aρ_0) indicates the average of A in the unperturbed equilibrium state ρ_0 and χ_AB(t-t') is known as the linear response function. 
    For the remainder of this work we will assume that the equilibrium density matrix ρ_0 commutes with H_0, and so is time independent. 
    The linear response function χ_AB(t-t') is given by the Kubo formula as
    
        χ_AB(t-t') = -iΘ(t-t')⟨[ A(t-t'),B(0)]⟩_0,
    
    where Θ(t-t') is the Heaviside step function, and the time evolution of operators is evaluated in the Heisenberg picture with respect to the unperturbed Hamiltonian H_0. 
    Assuming that f(t) goes to zero infinitely far in the past (such that the system can be viewed as starting in the unperturbed ground state at t→-∞), we can use the convolution theorem to take the Fourier transform of Eq. (<ref>) to find
    
        δ⟨ A⟩(ω) = ∫ dt e^iω tδ⟨ A⟩(t) = χ_AB(ω) f(ω),
    
    with 
    
        χ_AB(ω)    = lim_ϵ→ 0 -i∫_0^∞ dt e^i(ω + i ϵ)t⟨[ A(t),B(0)]⟩_0 
           ≡-i∫_0^∞ dt e^iω^+t⟨[ A(t),B(0)]⟩_0.
    
    In the second line of Eq. (<ref>) we have introduced the shorthand ω^+ = ω+iϵ, where it is understood that the limit ϵ→ 0 should be taken at the end of any calculation.
    
    We now review some textbook properties of χ_AB(ω), which can be found in, e.g. 
    Refs. <cit.>. 
    First, note that the Heaviside function in Eq. (<ref>) enforces causality: the response δ⟨ A⟩(t) depends only on the value of the perturbing field at earlier times t' < t. 
    Writing
    
        χ_AB(t-t') =  1/2π∫ dωχ_AB(ω) e^-iω(t-t'),
    
    this implies that χ_AB(ω) viewed as a function of complex ω_a+iω_b must be analytic in the upper half plane ω_b≥ 0. 
    To see this, we can imagine evaluating the inverse Fourier transform in Eq. (<ref>) via contour integration. 
    For (t-t')<0 we can close the contour in the upper half plane, since then the exponential factor is decaying for large |ω|. 
    This is guaranteed to give zero if there are no poles of χ_AB(ω) in the upper half plane. 
    We can verify this directly from Eq. (<ref>) by inserting a complete set of eigenstates {|n⟩} of H_0 (with energies E_n) into the average. 
    Taking the equilibrium density matrix to be given by Eq. (<ref>), 
    we have
    
        χ_AB(ω)    =  -i∫_0^∞ e^iω^+ t∑_nm p_n(⟨n|A|m⟩e^i(E_n-E_m)t⟨m|B|n⟩ -  ⟨n|B|m⟩⟨m|A|n⟩e^i(E_m-E_n)t) 
           =∑_nm p_n (⟨n|A|m⟩⟨m|B|n⟩/ω^+-(E_m-E_n) - ⟨n|B|m⟩⟨m|A|n⟩/ω^+-(E_n-E_m)).
    
    Eq. (<ref>) is known as the spectral (or Lehmann) representation of χ_AB(ω). 
    We see explicitly that the poles of χ_AB in the complex ω plane occur when ω = E_n-E_m-iϵ for every pair of energies E_n and E_m, and so are in the lower half plane as required by causality. 
    
    Further exploiting the analytic properties of χ_AB(ω) allows us to introduce a spectral density and derive Kramers-Kronig relations. 
    We will take a somewhat unconventional path to these results. 
    First, we derive an expression for Θ(t)e^iω^+ t, which appears in the integrand of Eq. (<ref>):
    
        Θ(t)e^iω^+ t = 1/2π i∫  dαe^iα t/α-ω^+.
    
    To verify this relation, we can evaluate the integral using contour integration in the complex α plane. 
    When t>0, we can choose to close the integration contour in the upper half complex α plane (so that the exponential factor decays at infinity), and so pick up the residue at the pole α=ω+iϵ. 
    On the other hand, when t<0, we must close the contour in the lower half plane, yielding zero since the integrand is analytic in the lower half plane. 
    Using this relation, we can write
    
        χ_AB(ω) = 1/π∫ dαχ_AB”(α)/α-ω^+,
    
    where χ_AB”(α) is known as the spectral density and given by 
    
        χ_AB”(α) = -1/2∫ dt e^iα t⟨[ A(t),B(0)]⟩_0.
    
    By inserting a complete set of states as we did in deriving the Lehmann representation Eq. (<ref>) of χ_AB, we can derive that the spectral density is given by
    
        χ_AB”(α) = π∑_nm⟨n|A|m⟩⟨m|B|n⟩δ(α+E_n-E_m) (p_m-p_n).
    
    We see that the spectral function is given by a weighted sum of δ-functions. 
    Note that when we consider response functions with A=B, then χ_AA” is real-valued. 
    This is not true in general, however.
    
    Next, we can use the Plemelj formula
    
        1/α-ω^+ = P1/α-ω + iπδ(α-ω)
    
    (where P denotes the Cauchy principal value) to derive several useful properties of χ_AB and χ_AB”. 
    Inserting Eq. (<ref>) into Eq. (<ref>) we find
    
        χ_AB(ω)    = 1/πP∫ dαχ_AB”(α)/α-ω + i χ_AB”(ω) 
           ≡χ_AB'(ω) + i χ_AB”(ω),
    
    where we have defined χ_AB'(ω) as the Hilbert transform of χ_AB”(ω). 
    When we consider response functions with A=B, we have from Eq. (<ref>) that χ_AA'(ω) is the real part of χ_AA(ω), and χ_AA”(ω) is the imaginary part of χ_AA(ω). 
    For more general response functions, this does not hold.
    
    By exploiting the analyticity of χ_AB(ω) in the upper half plane, we also have that
    
        0 = 1/2π i∫ dαχ_AB(α)/α-ω^-.
    
    Using the Plemelj formula Eq. (<ref>) on the right hand side and comparing with Eqs. (<ref>) and (<ref>) yields the Kramers-Kronig relations
    
        χ_AB'(ω)    = 1/πP∫ dαχ_AB”(α)/α-ω, 
        χ_AB”(ω)    = -1/πP∫ dαχ_AB'(α)/α-ω.
    
    
     §.§ Sum Rules
    
    
    Using Eq. (<ref>) we can derive an asymptotic expansion for χ_AB(ω→∞) in terms of moments of the spectral density χ_AB”(α). 
    We can obtain an asymptotic expansion by taking ω→∞ and Taylor expanding the denominator in Eq. (<ref>). 
    We find
    
        χ_AB(ω→∞)   ∼1/π∫ dα∑_n=0^∞-α^n/ω^n+1χ_AB”(α) 
           ∼∑_n=0^∞1/ω^n+1[-1/π∫ dαα^nχ_AB”(ω)]
           ∼∑_n=0^∞μ_AB^(n)/ω^n+1.
    
    From Eq. (<ref>) we see that the n-th moments μ_AB^(n) of χ_AB” determine the large-frequency asymptotic decay of the response function χ_AB(ω→∞). 
    Furthermore, since χ”_AB(α) is given in Eq. (<ref>) as the Fourier transform of a correlation function, we can integrate by parts to express μ_AB^(n) as an equal time average
    
        μ_AB^(n)   = -1/π∫  dαα^nχ”_AB(α) 
           =1/2π∫  dα   dt α^n e^iα t⟨[ A(t),B(0)]⟩_0
           =⟨[(i∂_t)^n A(0),B(0)]⟩_0
           =(-1)^n⟨[[H,[H,…[H,A]]]_n times,B]⟩_0.
    
    Eq. (<ref>) gives a family of sum rules, relating the n-th moment of the spectral function to the equal-time average of the commutator of the n-th time derivative of A with B. 
    Hence the sum rule also expresses the asymptotic decay of the response function χ_AB at large ω to the same equal time commutator, via Eq. (<ref>). 
    Importantly, Eqs. (<ref>) and (<ref>) are exact relations valid for any quantum system and any response function. 
    They hold for both noninteracting and strongly interacting electron systems, and represent one of the few exact constraints on measurable quantities.
    
    
    
     §.§ Application: Density-Density Response
    
    
    As an illustrative example, let us consider the density-density response function
    
        χ_ρρ(𝐪,ω) ≡χ_ρ_𝐪ρ_-𝐪(ω)
    
    that governs the response of the Fourier component ρ_𝐪 of the average density with wavevector 𝐪 to perturbations that couple to the particle density. 
    The corresponding spectral density is given by
    
        χ_ρρ”(𝐪,α) = -1/2∫ dt e^iα t⟨[ρ_𝐪(t),ρ_-𝐪(0)]⟩_0.
    
    We can now evaluate the moments μ_ρρ^(n)(𝐪,ω) of the spectral function χ_ρρ”. 
    We will focus on the lowest two cases, n=0 and n=1. 
    For n=0 we find
    
        μ_ρρ^(0)(𝐪,ω)    = -1/π∫  dαχ_ρρ”(𝐪,ω)
           = ⟨[ρ_𝐪,ρ_-𝐪]⟩_0 
           = 0,
    
    where we have used the fact that the density operators commute with each other at equal time, since they are functions only of the position operator (note that this relation fails to hold for the density operator projected into a set of low-energy bands). 
    Turning next to n=1, we have
    
        μ_ρρ^(1)(𝐪)    = -1/π∫  dααχ_ρρ”(𝐪,ω)
           = ⟨[i∂_tρ_𝐪,ρ_-𝐪]⟩_0.
    
    To simplify this further, we recall that charge conservation implies
    
        i∂_tρ_𝐪 = -[H,ρ_𝐪] = 𝐪·𝐣_𝐪,
    
    where 𝐣_𝐪 is the 𝐪-th Fourier component of the current density operator. 
    Inserting this into Eq. (<ref>) we find
    
        μ_ρρ^(1)(𝐪) = ⟨[𝐪·𝐣_𝐪,ρ_-𝐪]⟩_0.
    
    Eq. (<ref>) is the general form of the f-sum rule. 
    For nonrelativistic (n.r.) systems where the current density is proportional to the momentum density, the canonical commutation relations allow us to write the right hand side as 
    
        μ_ρρ^(1)(𝐪)→_n.r.|𝐪|^2n̅/m,
    
    where m is the particle mass and n̅ is the average ground state density n̅ = ⟨ρ_𝐪=0⟩_0.
    
    More generally, we can follow Ref. <cit.> to relate the average in Eq. (<ref>) to the diamagnetic current. 
    Concretely, note that for any operator 𝒪({𝐱̂_i},{𝐩̂_̂î}) 
    that depends on the position and momentum operators for each particle in the system (indexed by i), we can define an operator 𝒪_A = 𝒪({𝐱̂_i},{𝐩̂_i-𝐀(𝐱̂_i)}) 
    that is minimally coupled to a background electromagnetic vector potential 𝐀. 
    Writing the density operator as
    
        ρ_𝐪 = 1/V∑_i e^-i𝐪·𝐱_i,
    
    where V is the volume of the system, we have using the canonical commutation relations that
    
        [𝒪,ρ_𝐪] = . δ𝒪_A/δ A^μ_-𝐪|_𝐀→ 0q_μ.
    
    Applying Eq. (<ref>) to both the continuity equation Eq. (<ref>) and to the sum rule Eq. (<ref>) we find
    
        μ_ρρ^(1)(𝐪) = q_μ q_ν⟨.δ^2 H_A/δ A^μ_𝐪 A^ν_-𝐪|_𝐀=0⟩_0,
    
    where H_A is the Hamiltonian H_0 minimally coupled to the electromagnetic vector potential 𝐀. 
    Eq. (<ref>) shows that the first moment μ_ρρ^(1)(𝐪) of the spectral density χ_ρρ”(𝐪,ω) is proportional to the average of the diamagnetic current 
    
        .δ^2 H_A/δ A^μ_𝐪 A^ν_-𝐪|_𝐀=0.
    
    Combining Eqs. (<ref>), (<ref>), (<ref>), and (<ref>), we have that the high-frequency asymptotic expansion of the density-density response function is given by
    
        χ_ρρ(𝐪,ω→∞) ∼q_μ q_ν/ω^2⟨.δ^2 H_A/δ A^μ_𝐪 A^ν_-𝐪|_𝐀=0⟩_0.
    
    
    § DERIVATION OF THE SECOND-ORDER KUBO FORMULA
    
    
    In this section, we derive expressions for the linear and second order response functions for the average of an operator A in response to a perturbation
    
        H_1(t) = g(t)B.
    
    We will assume for now that the form of the operator A does not depend of g(t), so that we do not have any “contact” (diamagnetic) terms. 
    We assume that at large negative times t→ -∞ that the system is in a steady state of the unperturbed Hamiltonian H_0 with density matrix ρ_0 such that [H_0,ρ_0]=0. 
    We also assume that g(t) goes to zero as t→ -∞, such that we can write 
    
        g(t) = e^ϵ tf(t)
    
    with ϵ a positive infinitesimal that is taken to zero at the end of any calculation. 
    Eq. (<ref>) ensures that the perturbation f(t) turns on adiabatically. 
    
    We define the linear and second order response functions χ_AB(t-t') and χ^(2)_AB(t-t',t-t”) via a series expansion of ⟨ A⟩(t) in powers of f(t). 
    Concretely,
    
        ⟨ A⟩(t) = ⟨ A⟩_0 + ∫  dt' χ_AB(t-t')f(t') + ∫  dt'∫ dt”χ^(2)_AB(t-t',t-t”)f(t')f(t”)+….
    
    To evaluate the response functions, we can solve the equations of motion for the density matrix ρ(t) characterizing the state of the system in the presence of the perturbation H_1, subject to the initial conditions ρ(t→-∞)=ρ_0. 
    Writing
    
        ρ(t) = ρ_0 + ρ_1(t) +ρ_2(t)+…,
    
    where the subscript indicates the order of f on which each term depends, we can write the equations of motion for ρ(t) as
    
        ∂/∂ t   [ρ_0 + ρ_1(t) +ρ_2(t)+…] = -i[H,ρ] 
           =-i[H_0,ρ_1] -i[H_1,ρ_1] - i[H_0,ρ_2]+….
    
    Equating powers of f(t), we arrive at the following equations for ρ_1(t) and ρ_2(t):
    
        ∂ρ_1/∂ t   =-iH_0ρ_1-iH_1ρ_0,
        ∂ρ_2/∂ t   =-iH_0ρ_2-iH_1ρ_1.
    
    We can first solve Eq. (<ref>) by switching to the interaction picture, following Ref. <cit.>. 
    We find
    
        ρ_1(t)    = -i∫_-∞^t dt' e^i(t'-t)H_0H_1(t')e^-i(t'-t)H_0ρ_0
           =-i∫_-∞^t dt'B(t'-t)ρ_0f(t')e^ϵ t',
    
    where the time dependence in Eq. <ref> is evaluated using the unperturbed Hamiltonian H_0. 
    Inserting Eq. (<ref>) into Eq. (<ref>), we can solve for ρ_2(t) in the same way. 
    We find that
    
        ρ_2(t)    = -i∫_-∞^t dt' e^-iH_0te^it'H_0H_1(t')e^-it'H_0e^iH_0t'ρ_1(t')e^-iH_0t'e^iH_0t
           = -∫_-∞^t dt'∫_-∞^t' dt”B(t'-t)B(t”-t)ρ_0f(t')f(t”)e^ϵ(t'+t”).
    
    
    
    Using Eq. (<ref>) to write the average ⟨ A⟩(t) as
    
        ⟨ A⟩(t) = ⟨ A⟩_0 + tr(Aρ_1(t)) +tr(Aρ_2(t))+…,
    
    we can equate terms with Eq. (<ref>) order by order in f(t) to extract the response functions. 
    For the linear response function
    
        ∫  dt' χ_AB(t-t')f(t')    = tr(Aρ_1(t)) 
           = -i∫_-∞^t dt'tr(AB(t'-t)ρ_0)f(t')e^ϵ t'
           =-i∫_-∞^t dt'⟨A(0)B(t'-t)⟩_0 f(t')e^ϵ t'
           =-ie^ϵ t∫Θ(t-t')⟨A(t)B(t')⟩_0 f(t')e^-ϵ (t-t'),
    
    where we have used time-translation invariance to shift the time arguments within the average by t. 
    Eq. (<ref>) yields immediately Eq. (<ref>) of the main text.
    For the second-order response, we have
    
        ∫   dt' dt”χ^(2)_AB(t-t',t-t”)f(t')f(t”) =tr(Aρ_2(t)) 
           =  -∫_-∞^t dt'∫_-∞^t' dt”tr(AB(t'-t)B(t”-t)ρ_0)f(t')f(t”)e^ϵ(t'+t”)
           =-e^2ϵ t∫ dt'  dt”Θ(t-t')Θ(t'-t”)⟨A(t)B(t')B(t”)⟩_0e^-ϵ(2t-t'-t”)f(t')f(t”).
    
    We thus find that
    
        χ^(2)_AB(t-t',t-t”) = -Θ(t-t')Θ(t'-t”)⟨A(t)B(t')B(t”)⟩_0e^-ϵ(2t-t'-t”).
    
    As in linear response, it will be useful to Fourier transform Eq. (<ref>) to arrive at an expression for the frequency-dependent nonlinear response. 
    Concretely, we find
    
        1/2π∫  dt e^iω ttr(ρ_2(t)A)e^-2ϵ t   =1/2π∫  dt   dt' dt”χ^(2)_AB(t-t',t-t”)f(t')f(t”)e^iω t
           =1/2π∫  dt  dt' dt”  dω_1  dω_2 χ^(2)_AB(t-t',t-t”)f(ω_1)f(ω_2)e^iω te^-iω_1 t' e^-iω_2t”
           =∫  dω_1 dω_2 δ(ω-ω_1-ω_2)χ^(2)_AB(ω_1,ω_2)f(ω_1)f(ω_2),
    
    where we have defined
    
        χ^(2)_AB(ω_1,ω_2)    = ∫ dx∫ dy χ^(2)_AB(t_1,t_2) e^iω_1 t_1e^iω_2 t_2
           =-∫ dt_1∫ dt_2 Θ(t_1)Θ(t_2-t_1)⟨A(t_2)B(t_2-t_1)B(0)⟩_0e^iω_1^+ t_1e^iω_2^+ t_2.
    
    Note that because the operators A and B do not necessarily commute with themselves nor each other at different times, we cannot extend the domain of time integration by symmetrizing under ω_1↔ω_2 (in contrast to the case of time-ordered correlation functions).
    
    To make the causal structure of χ^(2)_AB(ω_1,ω_2) manifest, it is useful to perform an additional change of variables
    
        u   =t_1,
        
        v   =t_2-t_1
    
    in Eq. (<ref>) to find
    
        χ^(2)_AB(ω_1,ω_2)=-∫_0^∞du∫_0^∞dv ⟨A(u+v)B(v)B(0)⟩_0e^iω_12^+ ue^iω_2^+ v,
    
    which is the form of χ^(2)_AB(ω_1,ω_2) introduced in the main text Eq. (<ref>).
    
    To conclude, let us note that causality places analyticity constraints on χ^(2)_AB(ω_1,ω_2). 
    First, note from Eq. (<ref>) that χ^(2)_AB(t_1,t_2) vanishes for t_1<0 and for t_2<0 (in fact, it vanishes whenever t_2<t_1). 
    Writing
    
        χ^(2)_AB(t_1,t_2) = 1/4π^2∫ dω_1∫ dω_2 χ^(2)_AB(ω_1,ω_2)e^-i(ω_1 t_1+ω_2t_2)
    
    and evaluating the integrals by contour integration, we deduce that χ^(2)_AB(ω_1,ω_2) is analytic when both Im(ω_1)>0 and Im(ω_2)>0.
