A mathematical analysis of the adiabatic Dyson equation from time-dependent density functional theory

In this article, we analyze the Dyson equation for the density-density response function (DDRF) that plays a central role in linear response time-dependent density functional theory (LR-TDDFT). First, we present a functional analytic setting that allows for a unified treatment of the Dyson equation with general adiabatic approximations for discrete (finite and infinite) and continuum systems. In this setting, we derive a representation formula for the solution of the Dyson equation in terms of an operator version of the Casida matrix. While the Casida matrix is well-known in the physics literature, its general formulation as an (unbounded) operator in the N-body wavefunction space appears to be new. Moreover, we derive several consequences of the solution formula obtained here; in particular, we discuss the stability of the solution and characterize the maximal meromorphic extension of its Fourier transform. We then show that for adiabatic approximations satisfying a suitable compactness condition, the maximal domains of meromorphic continuation of the initial density-density response function and the solution of the Dyson equation are the same. The results derived here apply to widely used adiabatic approximations such as (but not limited to) the random phase approximation (RPA) and the adiabatic local density approximation (ALDA). In particular, these results show that neither of these approximations can shift the ionization threshold of the Kohn-Sham system.


Introduction 2 Introduction
Time-dependent density functional theory (TDDFT) is a formally exact theory to study the time evolution of a system of electrons; it has many applications in quantum chemistry, condensed-matter physics, and material science [7,11,35,36].Most of these applications lie within the perturbative regime, where linear response theory applies (LR-TDDFT).In this regime, one is no longer interested in the whole non-linear evolution of the single-particle density of the system but instead in the linear dynamical response of the density to a variation of the external potential.Stated differently, one is interested in the density-density response function (DDRF) of the system.The fundamental equation of LR-TDDFT is the celebrated Dyson equation that formally connects the DDRF of a given system of interest to the DDRF of an equivalent system of non-interacting electrons, the Kohn-Sham system.The equivalence is in the sense that both systems have the same ground state density.In shorthand notation, the Dyson equation reads where ⋆ denotes the convolution in time, χ is the DDRF of the system of interest, χ 0 is the DDRF of the Kohn-Sham system, and F Hxc is the linear operator whose Schwartz kernel is the Hartree plus exchange-correlation (Hxc-)kernel of TDDFT.In principle, the Hxc-operator depends on the ground state density of the system or, equivalently, on the Kohn-Sham ground state density.In practice, this density dependence is highly non-trivial, and the exact Hxcoperator is unknown; one then relies on approximations of this operator.
While several classes of approximations for the Hxc-operator were suggested in the physics literature (see [21,35] for an overview), the overwhelming majority of calculations are performed with adiabatic approximations.In the adiabatic approximation, the Hxc-operator acts instantaneously on the perturbing potential and can be seen as an operator acting on the spatial part of the perturbation.Within the adiabatic approximation, the Dyson equation becomes χ F (t) = χ 0 (t) + t 0 χ 0 (t − s)F χ F (s)ds, (1.1) where the adiabatic approximation F is now an operator from the (tangent) space of densities to the (tangent) space of potentials, χ 0 is again the DDRF of the Kohn-Sham system, and χ F is now an approximation of the true DDRF of the system of interest.For suitable choices of F , such approximations are observed to reproduce many response properties of large quantum systems accurately (see, e.g., [36]).
In the non-perturbative regime, TDDFT models were considered in various settings (see, e.g., [13,9,10,18,32,26]).By contrast, the mathematical literature on the linear response regime is scarce.Some works [4,5,6,34] have focused on the numerical aspects of extracting relevant properties (e.g., excitation energies, oscillator strengths, and absorption spectrum cross-section) of the solution χ F in the finite-dimensional case, i.e., after discretization of the underlying function space.However, a first step towards a rigorous understanding of the Dyson equation (1.1) in the infinite-dimensional (continuum) setting has been taken only recently in [14].There, the authors presented a mathematical framework for studying the Dyson equation within the Random Phase Approximation (RPA).
The current paper extends the framework presented in [14] to deal with more general adiabatic approximations.Moreover, the main contributions of this article can be summarized as follows: (1) We generalize the functional analytic setting presented in [14] to allow for a unified treatment of the Dyson equation with general adiabatic approximations for both discrete (finite and infinite) and continuum systems.Notably, this new setting allows us to study the celebrated adiabatic local density approximation (ALDA).
(2) We derive an explicit representation formula for the solution χ F in the case where χ 0 is the DDRF of a self-adjoint Hamiltonian.For this, the key result is a representation formula for the Fourier transform χ F in terms of an operator version of the Casida matrix.
(3) We derive and discuss several consequences of this representation formula.In particular, we characterize the maximal meromorphic extension of χ F and show that, for widely used adiabatic approximations such as the RPA and ALDA, the maximal domain of meromorphic continuation of χ F and χ 0 are the same.Physically, this means that these approximations are not able to shift the ionization threshold of the Kohn-Sham system.
Remark (Response theory terminology).In the physics literature, the name density-density response function usually refers to the Schwartz kernel of χ H (t).Here we refer instead to the operator-valued function t → χ H (t) as the density-density response function.Let us also remark that χ H is sometimes called the (linear) susceptibility [7] or the reducible (or irreducible) polarizability operator [20].

Main results
We now introduce some notation and discuss our main results.Throughout this article, H is a self-adjoint operator acting on the anti-symmetric N-fold tensor product of L 2 (Ω, dµ), where (Ω, µ) is a measure space.The specific measure space is not relevant to our results; in particular, µ can be the counting measure on some countable set Ω ⊂ R n (discrete systems), the Lebesgue measure on some open set Ω ⊂ R n (continuum systems), or a combination of both (continuum systems with internal spin).Moreover, we assume the following.
Assumption 1.The self-adjoint operator H : D(H) ⊂ H N → H N satisfies the following: (i) (Spectral gap) The ground state energy (ii) (Real Hamiltonian) H commutes with complex conjugation.
Since the ground state of H is non-degenerate, we can unambiguously define its ground state single-particle density (or simply density) as where Ψ 0 is the unique (up to phase) normalized ground state wave function of H.We then introduce the norms and and define the respective weighted L 2 spaces as where supp(ρ 0 ) denotes the support of ρ 0 .As usual, we identify the functions that coincide µ-almost everywhere.Let us also introduce the reduced Hamiltonian where

The density-density response function
The density-density response function (DDRF) of H is the operator-valued function where θ(t) is the Heaviside step function, E 0 is the ground state energy of H, and the operators B = B Ψ 02 and B * = B * Ψ 0 are defined as follows: The connection of this definition with the linear response of the system associated with H can be found in [14] and will not be discussed here.Instead, we shall limit ourselves to describing the main properties of χ H and presenting a characterization of the maximal meromorphic extension of its Fourier transform.
Let us start by recalling some results from [14].In [14], it is shown that the DDRF of typical Schrödinger operators on R 3N is a uniformly bounded and strongly continuous function with values in the set of bounded operators between suitable L p spaces.By adapting the proof in [14,Section 2] to the current setting, we can show that χ H is, in fact, uniformly bounded and strongly continuous on the space of bounded linear operators from L 2 ρ 0 to L 2 1 /ρ 0 , denoted here by B(L 2 ρ 0 , L 2 1 /ρ 0 ).Consequently, the Fourier transform of χ H , defined via is a tempered distribution with values on B(L 2 ρ 0 , L 2 1 /ρ 0 ).(See, e.g., [14,Proposition 2.8] for a derivation of (2.10)).In particular, χ H has an analytic extension to the upper half-plane.With some abuse of notation, we denote this extension also by χ H .
An immediate consequence of the spectral gap assumption on H is that χ H can be analytic extended to the larger set C\ σ(H # ) ∪ σ(−H # ) .In fact, it is shown in [14] that χ H can be meromorphic extended to the domain where Γ := inf σ ess (H # ) is the ionization threshold of H.However, this extension is not maximal in general.For instance, some cancellations can occur due to the conjugation of the resolvent of H # with the operator B. To precisely characterize the maximal meromorphic extension of χ H , let us define the single-particle excitation spectrum3 of H as where B ǫ (λ) ⊂ C denotes the ball centered at λ with radius ǫ and P H # (U) denotes the spectral projection of H # on some Borel subset U ⊂ C. We also define the discrete and essential parts of σ 1 (H # ) as The first result of this paper is then the following. where In particular, the set of poles of χ H is given by σ disc 1 (H # ) ∪ σ disc 1 (−H # ).Moreover, these poles are all simple, and their rank is given by Remark (Poles with infinite rank).In fact, any isolated point in σ ess (H # ) can also be seen as a pole of infinite rank of χ H .The reason for excluding such poles from the discrete spectrum is that compact perturbations of the operator H # may turn these poles into an accumulation point of poles, thereby making these singularities no longer isolated.
Remark (Dark excitations).The complementary spectrum σ(H # ) \ σ 1 (H # ) corresponds to the dark excitations of H, i.e., the part of the spectrum that can not be obtained by shining light on the system and measuring the absorption cross-section.

Well-posedness of the Dyson equation
Let us now turn to the solutions of the adiabatic Dyson equation (1.1).The first order of business is to agree on the underlying solution space.In LR-TDDFT, the goal of the Dyson equation is to approximate the DDRF of a system of interacting electrons via the DDRF of the equivalent non-interacting Kohn-Sham system.The equivalence is in the sense that the Hamiltonians of both the interacting and non-interacting systems have the same ground state density ρ 0 .In particular, the DDRF of both systems should lie on the space of strongly continuous maps from R + to B(L 2 ρ 0 , L 2 1 /ρ 0 ), denoted here by Therefore, it seems natural to study the well-posedness of the Dyson equation within this space.This choice is not unique, and we shall further motivate it later.Nevertheless, the Dyson equation is well-posed in this space under a compatible boundedness assumption on the adiabatic approximation of the Hxc-operator.
. Then, there exists a unique solution χ F of the Dyson equation (1.1) in the space C s R + ; B(L 2 ρ 0 , L 2 1 /ρ 0 ) .Moreover, the solution map The proof of Theorem 2.2 is a standard application of Banach's fixed point theorem.For the details, we refer the reader to [14, Section 3], where the same theorem in a different function space is proved.Although the proof is rather simple, we show that Theorem 2.2 guarantees the well-posedness of the Dyson equation for widely used adiabatic approximations of the Hxc-operator under general conditions on the ground state density ρ 0 .In addition, the bijectivity of the solution map implies that, for any F ∈ B(L 2 1 /ρ 0 , L 2 ρ 0 ), the DDRF of a Hamiltonian with ground state density ρ 0 can be obtained by solving the Dyson equation for a unique reference χ 0 .Of course, this does not guarantee that χ 0 is the DDRF of a non-interacting Hamiltonian, a common premise of LR-TDDFT.

The Dyson density-density response function
Throughout this section, we assume that χ F is the unique solution of the Dyson equation (1.1) with χ 0 = χ H for some H satisfying Assumption 1.Our goal is then to characterize these solutions and establish some of their fundamental properties.For this, the key ingredient is a representation formula for χ F based on an operator version of the Casida matrix.More precisely, let us formally define the Casida operator as Then under the assumption that F is symmetric (which is satisfied for physically relevant adiabatic approximations), we can apply the KLMN theorem (see Section 3) to properly define C as a semi-bounded self-adjoint operator on {Ψ 0 } ⊥ .Moreover, one can show (see Lemma 4.3) that defines a bounded operator on {Ψ 0 } ⊥ .The key result of this paper is that the Fourier transform of χ F is given by the conjugation of B with the operator above.Precisely, we have Theorem 2.3 (Solution in the frequency domain).Let F ∈ B(L 2 1 /ρ 0 , L 2 ρ 0 ) be symmetric, then the Fourier transform of χ F is well-defined for |Im(z)| > B * F B and satisfies (2.14) We can now derive several properties of the solution χ F .For starters, we can take the inverse Fourier transform of eq.(2.14) to obtain the following representation formula for χ F in the time domain.
Corollary 2.4 (Solution in the time domain).The solution χ F is given by the formula where sinc(s) = sin(s)/s is the analytic sinc function.
Remark.As the sinc function has a power series with only quadratic terms, the function sinc( √ s) defines an entire function.Moreover, this function is bounded on any half-line [α, ∞).In particular, sinc(t √ C) uniquely defines a bounded operator for any self-adjoint operator bounded from below.The fact that H The above representation formula highlights some important features of the solution χ F .For instance, note that we can decompose where C + and C − are the non-negative self-adjoint operators corresponding respectively to the positive real part (on (0, ∞)) and the non-negative part (on (−∞, 0]) of the spectrum of C. As these are mutually orthogonal operators, we have the decomposition .
Consequently, if 0 ∈ σ(C), then the positive part χ + F is stable in the sense that it is uniformly bounded in time, as expected from a DDRF of an isolated quantum system.The negative part, on the other hand, grows exponentially fast with time.Nevertheless, note that, since the operator C is bounded from below, the exponential growth of χ − F is bounded by When 0 ∈ σ(C), the solution may contain a linearly growing part, corresponding to the spectral projection of C on 0. In particular, if we gradually increase the strength of the adiabatic approximation by setting F (ǫ) = ǫF , the point ǫ 0 where the spectrum of the Casida operator reaches 0 corresponds to a phase transition of the system.The next corollary shows that the stability condition C > δ can be re-stated in terms of the simpler operator (2.16) Corollary 2.5 (Stability condition).Let M be the operator defined in (2.16), then we have In particular, the solution χ F is stable (in the sense described above) if and only if for some δ > 0.Moreover, χ F is a tempered distribution and χ F admits an analytic extension to the upper half-plane if and only if Since χ F is supposed to approximate the DDRF of another Hamiltonian (which is causal and bounded), the stability condition is expected to hold for physically relevant F .For the RPA, condition (2.17) follows from the positivity of F RPA (see eq. (2.21)) and the fact that inf σ(H # ) > 0. For the ALDA, however, we are not aware of a general argument to prove that (2.17) or (2.18) holds.

Remark (Quantitative stability). A quantitative version of Corollary 2.5 can be obtained via the min-max principle. Specifically, one can show that
where Remark (Stability in the finite-dimensional case).A finite-dimensional analog of M, and the associated stability condition, also appear in previous works where linear-response eigenvalue problems in finite dimensions are considered [33,22,4,5,6].More precisely, M is an operator version of the matrix M = A + B defined in the space of orbital pairs in [6].
Theorem 2.3 also allows us to characterize the maximal meromorphic extension of χ F .This characterization requires a spectral gap assumption on C and resembles the characterization of χ H given in Theorem 2.1.To state it precisely, let us define the single-particle spectrum of C as # P C (B ǫ (λ)) = 0 for any ǫ > 0 small}, where P C is the spectral projection of C. As before, we define the discrete part of σ 1 (C) as the set of isolated points with rank BH 1 /2 # P C ({λ}) < ∞, and the essential part as the complement of the discrete part.Then, we have the following characterization.
Corollary 2.6 (Maximal meromorphic extension of χ F ). Suppose that σ 1 (C) has a spectral gap on the non-negative part of the spectrum, i.e., [0, ∞) ⊂ σ 1 (C).Then the maximal meromorphic extension of χ F is given by where In particular, the set of poles of χ F is given by {z ∈ C : z 2 ∈ σ disc 1 (C)}.Remark (Poles rank and order).The non-zero poles of χ F are all simple and satsify On the other hand if 0 ∈ σ disc 1 (C), then χ F has a pole of second order at the origin.As a last consequence of Theorem 2.3, we show that under a compactness condition on F , the domain of maximal meromorphic continuation of χ F agrees with the domain of maximal meromorphic continuation of χ H .As the proof of this result is more involved than the proof of the previous corollaries, we promote it to a theorem.
The relevance of the above theorem is that for typical Schrödinger operators, the maximal domain of meromorphic continuation of χ H is related to its ionization threshold via eq.(2.11).Since the compactness condition (2.20) holds for standard adiabatic approximations (see Proposition 2.9 below), the above corollary implies that such approximations are not able to shift the ionization threshold of H (which is the Kohn-Sham Hamiltonian in applications).

Applications
We now discuss some applications of the previous results in the context of LR-TDDFT.Throughout this section, we work within the quantum chemistry set-up where Ω = R 3 and the underlying single-particle space is the classical Lebesgue space L 2 (R 3 ).
In this setting, the typical Hamiltonians of interest (e.g., the molecular Hamiltonian) are Schrödinger operators where ∆ is the Laplacian on R 3N and V is some real-valued function that acts by multiplication.Under general assumptions on V (e.g., V is in the Kato class of R 3N [31]), the ground state density of H is bounded whenever it exists.In particular, the following criterion applies to many situations encountered in practice.(The proof is a straightforward application of Hölder's inequality.) . The above criterion is easily verified for the following adiabatic approximations: • The random phase approximation (RPA).In the RPA, F is given by Thus from the Hardy-Littlewood-Sobolev (HLS) inequality, we conclude that F RPA ∈ B(L 2 1 /ρ 0 , L 2 ρ 0 ).(In fact, we just need ρ 0 ∈ L 3 2 (R 3 ) here.) • The Petersilka, Gossmann, and Gross approximation (PGG) [24].In the PGG approximation, the operator F is given by where γ H (r, r ′ ) is the ground state single-particle density matrix of the Hamiltonian associated to χ H . Hence, from the simple inequality |γ H (r, r ′ )| 2 ≤ ρ 0 (r)ρ 0 (r ′ ) and the HLS inequality, we also have • The adiabatic local density approximation (ALDA) [37,21,35].The ALDA is not a single approximation but rather a class of approximations.In the ALDA, the operator F is given by , is the exchange-correlation energy per particle of the homogeneous electron gas.While the exchange part is known and given by the correlation can only be approximated, which leads to different approximations of F ALDA ρ 0 . To see why such approximations also belong to B(L 2 1 /ρ 0 , L 2 ρ 0 ), let us take the parametrization of ε HEG c introduced by Perdew and Wang [23] as an example.The PW correlation approximation is where P = 1 or 3 4 , and A, α 1 , β 1 , β 2 , β 3 , β 4 > 0 are parameters chosen to reproduce the asymptotics expansions of ε HEG c in the low and high-density limits, and to fit data from quantum Monte Carlo simulations [12] in the intermediate regime.Thus from (2.22) and (2.23) (and some tedious calculations), one can check that Therefore, for any bounded ρ 0 .Other parametrizations of ε HEG c (ρ) also satisfy the above inequality as long as they reproduce (up to second derivatives) the asymptotic expansion of ε HEG c in the low-density limit.
The above list contains the most common adiabatic approximations used in practice and is not exhaustive.Note also that all adiabatic approximations mentioned above are symmetric.In particular, Proposition 2.8 guarantees that the solution formulas derived here apply to the Dyson equation with these approximations under the sole condition that the ground state density of the Kohn-Sham system is bounded.
Remark (Absorption spectrum).It turns out that the ground state density of typical quantum systems is not only bounded but also decays exponentially fast at infinity [1,2,31].In this case, the weighted density space L 2 ρ 0 contains functions that can grow exponentially fast at infinity.In particular, the polarizability tensor is a well-defined tempered distribution.Similarly, if the stability condition (2.18) holds, the polarizability tensor of the solution χ F also defines a tempered distribution, which can be shown to display peaks (Dirac's delta) at the poles of χ F (see [8,Proposition 17 and 44]).In applications, these peaks correspond to approximations of the electronic excitation spectrum of many-body quantum systems.Quantifying how good these approximations are is an important open problem.
As a final result, we present a simple sufficient criterion for the compactness property (2.20) that applies to the aforementioned adiabatic approximations.As a consequence, Theorem 2.7 shows that none of these adiabatic approximations are able to shift the ionization threshold of the Kohn-Sham system.Proposition 2.9 (Compactness criterion).Suppose that ρ 0 ∈ L 1 (R 3 ) ∩ L ∞ (R 3 ) and that the domain of H is continuously embedded in the classical Sobolev space 1 /ρ 0 and some δ > −1, (2.24) Remark (Optimality of (2.24)).The condition δ > −1 in Proposition 2.9 is optimal, as the following example shows.Let N = 1 and define the adiabatic approximation F ρ 0 as for some c ∈ R. Since N = 1, the ground state density is simply ρ 0 (r) = |Ψ 0 (r)| 2 and the operators B and B * reduce to Thus we have C = H 2 # + 2cH # , which implies that 2 , provided that c = 0 and σ ess 1 (H # ) is non-empty.Outline of the paper.We introduce some notation in the next paragraph.In the following section, we recall some well-known results about self-adjoint operators, their quadratic forms, and the associated Sobolev scale of spaces.These results are then used in Section 4 to prove the main results of this paper.In the appendix, we briefly discuss the optimality of weighted density spaces as the domain and co-domain of the DDRF.

Notation
We denote the set of non-negative real numbers by R + .For A and B scalar quantities, A B means that there is an irrelevant positive constant C such that |A| ≤ C|B|.Occasionally, we also use A ǫ B to indicate the dependence of the implicit constant on the additional parameter ǫ.
Let F be a Banach space, then we denote its norm by • F , or simply by • if the space is clear from the context.The set of linear continuous operators from F to another Banach space G is denoted by B(F, G); the set of compact operators is denoted by B ∞ (F, G).The operator norm is denoted by T F,G or simply by T if the Banach spaces are clear from the context.The kernel and the range of T are denoted, respectively, by ker T ⊂ F and ran T ⊂ G.We also use rank T = dim ran T for the rank of T .The anti-dual of a Banach space F , i.e., the space of antilinear continuous functions from F to C endowed with the operator norm is denoted by F ⋆ .For the Fourier transform of a function f : R → F , we use the physics convention (2.25) For any Hilbert space H, we adopt the convention that the inner-product •, • H is antilinear in the first variable and linear in the second.For 1 ≤ p ≤ ∞, L p (R n ) denotes the standard L p spaces with respect to the Lebesgue measure on R n .We also use L p (R n ) + L q (R n ) and L p (R n ) ∩ L q (R n ) for the Banach spaces of Lebesgue-measurable functions with the norms If K j = 0 for some j ≥ 1, we say that z 0 is a pole of K.In addition, if K j = 0 for j ≥ 2, then we say that z 0 is a simple pole of K and define its rank as Moreover, we say that D is the maximal domain of K if there exists no meromorphic extension of K to a strictly larger connected domain.

Mathematical background
In this section we briefly recall some well-known facts about the scale of Sobolev spaces associated to self-adjoint operators and their quadratic forms.We also use this short recap to set-up some additional notation that will be used during the proofs from Section 4. The material presented here can be found in standard references such as [27,28,29].

Sobolev scale of spaces
Throughout this section, we let A : D(A) ⊂ H → H be some self-adjoint operator on a Hilbert space H satisfying the inequality Ψ, Aψ ≥ Ψ 2 for any Ψ ∈ D(A).(In the case A ≥ 1 − c for some c > 0, we replace A by A + c everywhere in the discussion below.)Then by the spectral theorem, there exists a measure space (X, ν), an unitary map U : H → L 2 (X, dν), and a real-valued ν-measurable function a : and A acts on L 2 (X, dν) by multiplication by a, i.e., UAU * f (x) = a(x)f (x).
We then define the Sobolev spaces induced by A as follows.
Definition (Sobolev spaces).For s ≥ 0, the Sobolev space of order s is the set Moreover, the negative Sobolev space of order −s is defined as the anti-dual space of H s (A), i.e., the set antilinear and continuous} endowed with the operator norm.
By the Riesz representation theorem, we can isometrically identify H with its anti-dual H ⋆ via the Riesz map In this way, we have a natural chain of dense inclusions From (3.1), the operator A s restricted to H m (A) defines an isometric isomorphism between H m (A) and H m−s (A) for any , 0 ≤ s ≤ m.Consequently, the adjoint map induces an isometric isomorphism from H s−m (A) to H −m (A).In particular, by the chain of inclusions in (3.2) and using the commutation relations the operator A s can be uniquely extended to a continuous4 operator on the whole Sobolev scale of spaces Furthermore, the chain of inclusions in (3.2) also allow us to naturally define the operator for any s ∈ R and µ ∈ C. A useful consequence of the spectral theorem is that the spectrum of A is independent of the Sobolev scale used.More precisely, we have Lemma 3.1 (Inverse on Sobolev scale).Let µ ∈ C, then the operator µ − A is invertible in B(H s (A), H s−2 (A)) for some s ∈ R if and only if it is invertible for every s ∈ R.
Proof.First, we can use the Riesz representation theorem on L 2 (X, dν) to extend the unitary map (given by the spectral theorem) U : H → L 2 (X, dν) to U : H −s (A) → L 2 (X, a −s dν) for any s ≥ 0.More precisely, UT ∈ L 2 (X, a −s dν) is the unique function satisfying Hence, the operator U(µ − A)U * acts on L(X, a −s dν) as pointwise multiplication by µ − a(x) for any s ∈ R. One can now check that this operator is invertible if and only if |µ − a(x)| > δ µ-a.e. for some δ > 0, which is equivalent to µ − A being invertible on H.
Let us conclude this section with a distributional version of Stone's formula (cf.[27, Theorem VII.13]) that will be useful to establish the maximality of the meromorphic extensions from Theorem 2.1 and Corollary 2.6.

Lemma 3.2 (Stone's formula).
Let A be a semi-bounded self-adjoint operator and R A (z) = (z − A) −1 denote the resolvent of A. Then for any f ∈ C ∞ c (R) we have where P A λ is the spectral projection-valued measure of A and the convergence is in the operator norm on B(H s (A), H s+2 (A)) for any s ∈ R.
Proof.Since (µ ± iη − λ) −1 is uniformly bounded in λ ∈ R for η > 0 fixed, from Fubini's theorem we have where p η (µ) = η µ 2 +η 2 is the Poisson kernel.As g(A) commutes with A for any continuous function g, it is enough to show that This now follows from the continuity of the spectral calculus and the estimate which can be shown by using the Fourier transform of the Poisson kernel p η (ω) = πe −η|ω| .

Quadratic forms
Let us now introduce the quadratic form associated to a semi-bounded operator A and present the KLMN theorem [28, Theorem X.17].For a proof, consult [28].
Definition (Quadratic form).For a semi-bounded self-adjoint operator A satisfying A ≥ 1−c for some c ∈ R, the associated quadratic form is the sesquilinear map q A : H 1 (A)×H 1 (A) → C defined as The Sobolev space of order 1 is also called the form domain of A.
Note that, by definition, Next, let β : H 1 (A) × H 1 (A) → C be another symmetric sesquilinear form.Then we say that β is relatively bounded with respect to q A if there exists some 0 < a < 1 and b > 0 such that With this definition, we can now state the KLMN theorem, which is essentially a quadratic form version of the celebrated Kato-Rellich theorem.
Lemma 3.3 (KLMN theorem [28]).Let A be a self-adjoint operator satisfying A ≥ 1 − c and β be a symmetric sesquilinear form on H 1 (A) satisfying (3.3).Then, there exists an unique self-adjoint operator B with the same form domain as A and satisfying The proof of the KLMN theorem consists in using inequality (3.3) to show that the norms q B + α •, • and (A + c) 1 2 • are equivalent and exploit the one-to-one relation between semibounded closed quadratic forms and semi-bounded self-adjoint operatos.In particular, using interpolation theory one can show that for any Φ ∈ H s (A), (3.5) provided that −1 ≤ s ≤ 1. (The identity above means that these are the same subsets of H for s ≥ 0, and not simply isometric.)

Reducing subspaces and block decomposition
We now recall the definition of reducing subspaces for unbounded self-adjoint operators and the associated block decomposition.For the proofs of the results presented next, we refer to [30,Section 1.4].As before, we assume that A : In particular, χ H : R → B(L 2 ρ 0 , L 2 1 /ρ 0 ) is bounded and strongly continuous.We can now complete the proof of Theorem 2.1 Proof of Theorem 2.1.Since χ H is bounded and strongly continuous, its Fourier transform is well-defined as a tempered distribution.Moreover, since χ H is causal (i.e., χ H (t) = 0 for t ≤ 0), the Fourier transform is analytic on the upper half-plane {z ∈ C : Im(z) > 0}.From the spectral Theorem and straightforward computation (see [14,Section 2]), this analytic continuation is given by Now note that, since the single-particle excitation spectrum σ 1 (H # ) is closed (thus Borel measurable), we can decompose the resolvent of H # in two terms: From the definition of σ 1 (H # ), the second term vanishes when multiplied by B on the left.In particular, the spectral theorem yields Next, observe that the spectral gap assumption on H implies that the set {z ∈ C : ±z ∈ σ 1 (H # )} is open and connected.Thus the right-hand side of (4.2) defines the unique analytic extension of χ H to this set.Moreover, for any isolated point λ 0 ∈ σ 1 (H # ), we have where the second term is analytic around λ 0 .In particular, χ H is meromorphic on From the identity rank S = rank SS * , we also obtain the rank equality To conclude the proof, we need to show that D is the maximal domain of meromorphic continuation of χ H .So suppose that χ H is analytic around some λ 0 ∈ R, then it is enough to show that λ 0 ∈ σ 1 (H # ).Moreover, since χ H (z) = χ H (−z), we can assume (without loss of generality) that λ 0 ∈ R + .Then, define As the right-hand side is continuous close to λ 0 , Stone's formula (see Lemma 3.2) yields 0 = lim with support in B ǫ (λ 0 ) for ǫ > 0 small enough.Choosing a sequence f n converging monotonically to the indicator function on B ǫ (λ 0 ) ∩ R and using the strong convergence property of the spectral calculus (see [27,Theorem VIII.5.(d)]), we find that BP H # (B ǫ (λ 0 ))B * = 0. Therefore, λ 0 ∈ σ 1 (H # ), which completes the proof.
Remark (Regularity of χ H along the continuous spectrum).If H has compact resolvent (e.g., a Schrödinger operator with a trapping potential), then the maximal domain of meromorphic continuation is the whole complex plane.However, for typical Hamiltonians in electronic structure theory (e.g., the atomic and molecular Hamiltonians), the spectrum is divided into discrete and continuous parts [29].In some exceptional cases, and for suitable f and g, the regularity of the map ω → g, χ H (ω)f along the continuous spectrum can be rigorously studied (see [15]) via the celebrated limiting absorption principle [3,16].

Proof of Theorem 2.3
We split the proof of Theorem 2.3 into two parts.In the first part, we properly define the Casida operator and relate its resolvent to the inverse of the operator In the second part, we show that χ H is given by the conjugation of B with C(z) −1 via the convolution property of the Fourier transform and a well-known resolvent identity.

The Casida Operator
To properly define the Casida operator, we consider the quadratic form β : Recall that we assumed F ∈ B(L 2 1 /ρ 0 , L 2 ρ 0 ) to be symmetric, and where Ψ 0 and E 0 are the ground state and ground state energy of H. Therefore, we can use the KLMN theorem to prove the following proposition.
Proof.Note that q H 2 # = H•, H• is the quadratic form of the positive operator Consequently, to apply the KLMN theorem, we just need to check that β is relatively bounded with respect to q H 2 # .For this, we can use Cauchy-Schwarz and Young's inequality to obtain for any a > 0. To prove eq. ( 4.4) we can use eq.(4.5) and the simple estimate The next step is to relate the Casida operator to the operator Since H # is self-adjoint and the operators H −1 # and B * F B are both bounded and symmetric, the operator C(z) is closed and normal.Formally, the operator C is given by # .We thus expect that # in an appropriate sense.Rigorously clarifying this statement is the goal of the next lemma.
Lemma 4.3 (The resolvent of C).Let C be the Casida operator defined according to Proposition 4.2 and C(z) be the operator defined above.Then z 2 − C is invertible if and only if C(z) is invertible.In this case, the operator Proof.Let H 1 (C) be the first Sobolev space associated with C. Then by Lemma 3.1, a point z 2 is in the resolvent set of C if and only if the extension z 2 −C : 2), this is equivalent to the map # Φ being bijective (by the closed graph theorem).However, the right-hand side of the above is equal to In particular, we have Thus since H To show (4.7), we note that if either C(z) or z 2 − C is invertible, then from (4.8) we obtain Eq. (4.7) now follows by multiplying eq. ( # on the left and on the right.

Proof of Theorem 2.3
We are now in position to prove Theorem 2.3.For this, we shall use the following well-known resolvent identity.
Lemma 4.4 (First resolvent identity).Let C(z) be the operator defined in (4.6).Then, if the operators C(z) and Proof of Theorem 2.3.First, note that from the Dyson equation (1.1), the simple estimate sup t∈R χ H (t) ≤ B B * = B 2 , and Gronwall's inequality we have In particular, the Fourier transform is well-defined and analytic for Im(z) > B 2 F .In this case, by the convolution property of the Fourier transform, we have which is the frequency version of the Dyson equation.
The idea now is to find a representation formula for the inverse of the dielectric operator where the second equality comes from the expression For this, we can now use the resolvent identities from Lemma 4.4.Precisely, let C(z) be defined as in eq. ( 4 is invertible for any z with z 2 ∈ (0, ∞).We now claim that ε(z) −1 = 1 + 2BC(z) −1 B * F for any z with Im(z) large.
Indeed, from (4.10) we have On the other hand, eq.(4.11) implies that ε(z)(1 + 2BC(z) −1 B * ) = 1 as well, which proves our claim.To conclude, we note that from the (frequency) Dyson equation (4.12) we have Thus Theorem 2.3 follows from the equation above and Lemma 4.3.

Proof of Corollaries 2.4 to 2.6
We now present the proofs of the corollaries from Theorem 2.3.We start with the proof of Corollary 2.4.
Proof of Corollary 2.4.First, note that the power series of the sinc function is composed only of even terms.Thus the function defines an analytic function in the whole complex plane.Moreover, this function is bounded on any half-line {λ ∈ R : λ ≥ −c} for c > 0. In particular, the operator f (t 2 C) defined via the spectral theorem belongs to B({Ψ 0 } ⊥ ) for any t ∈ R. In fact, by taking any function g : R → [1, ∞) satisfying g(λ) ∼ √ λ for λ big we can re-write where h(t, λ) = g(λ)f (t 2 λ) is bounded on σ(C).So by recalling that H s (C) = H 2s (H # ) for any −1 ≤ s ≤ 1 (because H 1 (C) = D(H # ) with equivalence of norms), eq.(4.14) implies that Therefore, the operator-valued map defines a strongly continuous family of operators in B(L 2 ρ 0 , L 2 1 /ρ 0 ).To conclude the proof, we can use the identity to show that the Fourier transform of χ F (t)e −ηt and χ(t)e −iηt coincide for η large enough, and therefore, χ(t) = χ F (t) for every t ∈ R. (Note that equality holds everywhere because both functions are strongly continuous.) Next, let us turn to the proof of the stability criterion.
Proof of Corollary 2.5.First, note that where C(z) is the operator defined in (4.6).In particular, by Lemma 4.3 we have 0 ∈ σ(C) if and only if 0 ∈ σ(M).Moreover, from the proof of Lemma 4.3 we know that Since a self-adjoint operator is non-negative if and only if its quadratic form is non-negative, the above identity (and the fact that H Lastly, we present the proof of Corollary 2.6.
Proof of Corollary 2.6.Let us define Then by multiplying and dividing by g(C) 1 2 as we did in eq.(4.14), we see that χ(z) is bounded on B(L 2 ρ 0 , L 2 1 /ρ 0 ) for any z 2 ∈ σ 1 (C).Moreover, χ(z) is meromorphic on and its poles are located at {z ∈ C : z 2 ∈ σ disc 1 (C)}.From the spectral gap assumption on σ 1 (C), the open set D F is connected.Thus since χ F (z) = χ(z) for Im(z) big enough, the function χ is the unique meromorphic extension of χ F to D F .
To show that this extension is maximal, we can now use Stone's formula as we did in the proof of Theorem 2.1.Precisely, let µ 0 ∈ C with µ 2 0 ∈ R and suppose that χ F can be analytic extended to a neighborhood U of µ 0 .Then it suffices to show that µ 2 0 ∈ σ 1 (C).In the case µ 2 0 > 0, we can assume that U does not intersect the imaginary axis and define where we choose the branch of the square root such that lim η→0 α ± (η, µ) = 0. Note that α ± is continuous in a neighborhood of η = 0 and µ ∈ U. Thus from the continuity of χ F on U and Stone's formula in Lemma 3.2 we have 0 = lim As in the proof of Theorem 2.1, we can now use the strong-convergence property of the functional spectral calculus to conclude that µ 2 0 ∈ σ 1 (C).
In the case µ 2 0 < 0, we can choose the neighborhood U so that U does not intersect the real axis.By defining α ± (µ, η) via (4.17) with the opposite branch of the square root, and using similar arguments, one can prove that µ 2 0 ∈ σ 1 (C).Finally, if µ 0 = 0, then µ 0 is at most an isolated point in σ 1 (C) by the preceding arguments.In this case, however, µ 0 would be a pole of χ F , which is not possible because we assumed that χ F is bounded around µ 0 .Therefore, µ 2 0 ∈ σ 1 (C), which concludes the proof.

Proof of Theorem 2.7
For the proof of Theorem 2.7, it is more convenient to work with the single-particle excitation spectrum of H 2 # than of H # , where the former is defined as More precisely, note that from the definitions of D and D F (see (2.13) and (2.19)), the proof of Theorem 2.7 is done if we show that where the essential part of σ 1 (H 2 # ) is defined in the same way as for H # .For this, the first step is the following lemma.Then V H and V C are reducing subspaces for both H 2 # and C and we have

.20)
Proof.From the spectral theorem, V H and V C are clearly reducing subspaces for H 2 # and C, respectively.Let us then show that V H is reducing for C. First, we claim that To prove this claim, first note that there exists a sequence {(λ j , ǫ j )} j∈N ⊂ R × (0, 1] such that BP H 2 # (B ǫ j (λ j )) = 0 for any j and lim m→∞ m j=1 in the strong sense.As a consequence, for any Ψ, Φ ∈ D(H # ).From the first identity in (4.23), we see that (4.21) holds.From the last identity in (4.23), we find that P V ⊥ H maps D(C) to itself and V H and V ⊥ H are invariant subspaces for C. We thus conclude from Lemma 3.5 that V H is a reducing subspace for C.
Finally, to prove that V C is reducing for H 2 # and that the second inclusion in (4.20) holds, we can reverse the roles of C and H 2 # .More precisely, we let α > 0 be big enough (e.g., α > B * F B 2 + 1 will do) and note that where . So repeating the same steps from before with the roles of C and H 2 # exchanged, the result follows.Next, we use the compactness assumption from Theorem 2.7 to show that the essential spectrum of C| V H ,respectively, C| V C is equal to the essential spectrum of as a bounded operator from H −1 (H # ) to H 3 (H # ) for any µ > 0 big enough.In particular, by recalling the chain of inclusions in eq.(3.2), identity (4.26) holds in B({Ψ 0 } ⊥ ).Furthermore, from Lemma 4.5 we have the block decomposition Since B * F B ∈ B ∞ (D(H # ), {Ψ 0 } ⊥ ), the second term in (4.26) is compact in B({Ψ 0 } ⊥ ).Therefore, from Weyl's criterion, we conclude that Moreover, a similar argument shows that Eq. (4.25) now follows from the relations σ(f (A)) = f (σ(A)) and ker f (λ) − f (A) = ker λ − A (for injective f ), which is a well-known corollary of the spectral theorem.
We can now complete the proof of Theorem 2.7.
Proof of Theorem 2.7.Since σ 1 is closed, from the spectral theorem and the definitions of V C and V H we have # P C| V H ({λ}) ≤ rank P C| V H ({λ}) < ∞, which shows that λ ∈ σ ess 1 (C).The second inclusion in (4.27) follows from the same arguments with the roles of H 2 # and C interchanged.

Proof of Proposition 2.9
For the proof of Proposition 2.9, we shall use the following version of the Rellich-Kondrachov (aka compact Sobolev embedding) theorem.From the classical Sobolev embedding, the norms Φ j − Φ k L p are uniformly bounded (in j, k) for any 2 < p ≤ 2n n−2 .So choosing M > 0 arbitrarily large, the second term can be made arbitrarily small, which completes the proof.
Proof of Proposition 2.9.For the RPA, we have B * F RPA B = P Ψ ⊥ 0 T , where T is an integral operator with integral kernel given by T (r 1 , .., r N , r ′ 1 , ..., r ′ N ) = N .
In particular, from the assumption that D(H # ) is continuously embedded in H 1 (∆), it is enough to show that for any bounded sequence {Φ j } j∈N in H 1 (∆), there exists a subsequence that satisfies (after re-labeling the indices) To find such a subsequence, we apply the Rellich-Kondrachov theorem (Lemma 4.7).Precisely, let {Φ j } be the subsequence satisfying (4.28) and define I ǫ := {r ∈ R 3 : ρ 0 (r) > ǫ}.
Then by the Cauchy-Schwarz inequality (as in (4.1)) and assumption (2.24) we obtain Ψ 0 (r, r) Φ j (r, r) − Φ k (r, r) dr where B M is the ball of radius M centered at the origin.From the definition of I ǫ , the first term is bounded by ǫ 2δ+2 .Moreover, the R 3 -measure of I ǫ is finite (because ρ 0 ∈ L 1 ) and Lastly, we recall the definition of an operator-valued meromorphic function [16, App.C].Definition 2.10 (Meromorphic operator-valued function).Let D ⊂ C be open, then we say that K : D → B(F, G) is a meromorphic operator-valued function if for any z 0 ∈ D there exists (i) a neighborhood U z 0 ⊂ D of z 0 , (ii) finitely many finite rank operators {K j } j≤M ⊂ B(F, G), and (iii) a holomorphic function K 0