Hartree-Fock and Kohn-Sham time-dependent response theory in a second-quantization atomic-orbital formalism suitable for linear scaling

Hartree-Fock and Kohn-Sham time-dependent response theory in a second-quantization atomic-orbital formalism suitable for linear scaling

Thomas Kjærgaard,^1,a兲Poul Jørgensen,¹Jeppe Olsen,¹Sonia Coriani,^2,3and Trygve Helgaker³

1Lundbeck Foundation Center for Theoretical Chemistry, Department of Chemistry, Aarhus University, Langelandsgade 140, DK-8000 Århus C, Denmark

2Dipartimento di Scienze Chimiche, Università degli Studi di Trieste, via L. Giorgieri 1, I-34127 Trieste, Italy

3Centre for Theoretical and Computational Chemistry, Department of Chemistry, University of Oslo, P.O. Box 1033, N-1315 Blindern, Norway

共Received 14 May 2008; accepted 27 June 2008; published online 5 August 2008兲

We present a second-quantization based atomic-orbital method for the computation of time-dependent response functions within Hartree-Fock and Kohn-Sham density-functional theories. The method is suited for linear scaling. Illustrative results are presented for excitation energies, one- and two-photon transition moments, polarizabilities, and hyperpolarizabilities for hexagonal BN sheets with up to 180 atoms. ©2008 American Institute of Physics.

关DOI:10.1063/1.2961039兴

I. INTRODUCTION

Response-function theory describes how an observable responds when a molecular system is subjected to an external field. If the field is time dependent and oscillates with a given frequency, it causes the observable to oscillate and to become frequency dependent. The response of the observ- able may be expanded in powers of the field strength: the linear response of the system is determined by the linear response function, the quadratic response by the quadratic response function, and so on.¹

Molecular response properties may be calculated from the response functions by specifying operators for the ob- servable and for the applied fields. For example, the frequency-dependent polarizability and hyperpolarizability may be evaluated from the linear and quadratic response functions, respectively, by specifying the electric dipole- moment operator for both the observable and the oscillating electric fields 共lasers兲. From the poles and residues of the response functions, additional molecular properties can be obtained, including excitation energies to electronically ex- cited states, strength parameters for共multiphoton兲transitions to these states, and excited-state properties.² In view of the wealth of molecular information afforded by the response functions, it is important that they are evaluated efficiently.

This is further substantiated by the fact that in density- functional theory—no doubt the method of choice nowadays for large systems—information about electronically excited states is only available through the response functions as these states cannot be obtained directly by an energy optimization.

In this paper, we describe how molecular response func- tions may be derived for Kohn-Sham density-functional theory³ using an atomic-orbital 共AO兲 basis representation and the formalism of second quantization. In the derivation, we focus first on the response function for a Hartree-Fock state and then introduce the generalizations that are required for Kohn-Sham theory.

Time-dependent perturbation theory for Hartree-Fock states has a long history. The theory of first-order molecular properties was first given by Ball and McLacland, using the time-dependent Hartree-Fock approximation.^4,5Second- and higher-order molecular properties in Hartree-Fock theory were subsequently evaluated.^6,7Response-function theory as it is used today—where response functions are subjected to a pole and residue analysis to determine molecular properties for both ground and excited states and for transitions be- tween these states—was developed in the 1980s共Ref.2兲and has since been used to calculate a large variety of molecular properties. The quasienergy approach that was developed in the 1990s 共Refs. 8and9兲has tied response-function theory closely to energy-derivative techniques. Through the quasienergy approach, response functions may be straightfor- wardly derived for both variational and nonvariational wave functions and used to calculate molecular properties for both ground and excited states, as well as transition properties between these states.^2,9

The linear response eigenvalue equations exhibit, in gen- eral, a paired structure with pairs of eigenvalues with identi- cal magnitude and opposite signs. Furthermore, it is suffi- cient to calculate the eigenvector for one eigenvalue of such a pair, as the eigenvector for the other eigenvalue may be obtained by a simple rearrangement.² A method has previ-

a兲Electronic mail: [email protected].

0021-9606/2008/129共5兲/054106/23/$23.00 129, 054106-1 © 2008 American Institute of Physics

ously been developed to exploit this pairing to obtain a fast and reliable iterative scheme for solving linear response equations.¹⁰

Using an AO parametrization, Hartree-Fock and Kohn- Sham response-function theory was described by Larsen et al.¹¹in the language of first quantization. Here we present the second-quantization analog of the derivation by Larsen et al.¹¹ correcting at the same time some inconsistencies in their derivation concerning the pairing properties of the gen- eralized Hessian and metric matrices. The derivation of Larsenet al.¹¹was recently used by Corianiet al.¹²to obtain a linear-scaling algorithm for solving the response eigen- value and linear equations and for evaluating second-order molecular properties. We extend the linear-scaling develop- ment of Ref. 12 to the description of molecular properties that may be obtained from the quadratic response function and its residues. We furthermore analyze the involved matri- ces and eigenvectors in detail to show how the eigenvector pairing may be used to obtain effective iterative schemes in the AO basis.

A linearly scaling time-independent response theory has been presented by Niklasson, and co-workers^13,14 within a purification framework, and recently, Kussmann and Ochsen- feld have also reported time-dependent Hartree-Fock and Kohn-Sham calculations of the frequency-dependent polariz- ability and hyper-polarizability using a linear-scaling framework.¹⁵

The theoretical foundation for time-dependent perturba- tion theory as a means to obtain information about electroni- cally excited states in Kohn-Sham density-functional theory has been discussed by several authors^16–18and a large variety of molecular properties calculated—see, for example, Refs.

19–23. Accepting the standard approach, only minor modifi- cations must be made to generalize time-dependent Hartree- Fock theory to Kohn-Sham theory. Using the response- function approach, the generalization has been described in detail for properties that relate to linear, quadratic, and cubic response functions in Refs.24and25. We describe here the modifications that are required in an AO-based parameteriza- tion.

In the next section, we summarize those aspects of sec- ond quantization that are relevant for deriving response func- tions within a second-quantization AO-based formalism. In Sec. III, we identify response functions with terms in the expansion of the expectation value of a general observable.

Next, in Sec. IV, we derive equations for the first- and second-order corrections to the wave function. In Sec. V, we discuss the linear and quadratic response functions, their residues, and their poles—in particular, we examine the pair- ing of the excitation energies obtained from the linear re- sponse function. In Sec. VI, the expressions required for a linear-scaling implementation are derived. In Sec. VIII, the Hartree-Fock theory is generalized to Kohn-Sham theory. Fi- nally, features of the theory are illustrated by calculations on boron nitride.

II. SECOND-QUANTIZATION-BASED AO THEORY A. The density matrix

Consider a set of nonorthogonal atomic spin orbitals␹␮

with the Hermitian metric S. The creation and annihilation operators of the AOs fulfill the anticommutation relations²⁶

关a_␮^†,a_␯^†兴+= 0, 共1兲

关a_␮,a_␯兴+= 0, 共2兲

关a_␮^†,a_␯兴+=S_␯␮. 共3兲

For a single-determinant state 兩0典, the expectation values of the creation and annihilation operators are given by

⌬_␮␯=具0兩a_␮^†a_␯兩0典, 共4兲 where

兩0典=a_i^†a^†_j¯a_l^†兩vac典 共5兲

anda_i^†a_j^†¯a_l^†refer to the set of orthonormal molecular spin orbitals that are occupied in兩0典. Roman letters are used here for the orthonormal molecular spin orbitals and Greek letters for their atomic counterparts. The orthonormal spin orbitals may be expanded in the AO basis as

a_i^†=

兺

The elements of the density matrix D are in the AO repre- sentation given as

D_␮␯=

兺

i=1 n

C_␮iC_␯^*_i, 共7兲

where the summation is over all occupied spin orbitals and where the expansion coefficients are normalized over the metric:

C^†SC=I. 共8兲

The density matrixDis related to⌬by the transformation

⌬=S^TD^TS^T. 共9兲 To prove this transformation, we use

关a_p^†,a_␯兴+=

兺

关a_␮^†,a_p兴+=

兺

to show that

具vac兩al¯ajaia_␮^†a_␯a_i^†a^†_j¯a_l^†兩vac典

=具vac兩al¯aja_␮^†a_␯a^†_j¯a_l^†兩vac典+

兺

共11兲

Repeated use of Eq. 共11兲shows that Eq. 共9兲is valid. Thus, the AO density-matrix element D_␮␯ is only identical to the matrix element ⌬_␮␯in an orthonormal basis.

Turning our attention to the two-electron case, we show in Appendix A that the two-electron expectation value

⌫_␮␯␳␴=具0兩a_␮^†a_␳^†a_␴a_␯兩0典 共12兲 decouples into products of expectation values of one-electron operators,

⌫_␮␯␳␴=⌬_␮␯⌬_␳␴−⌬_␮␴⌬_␳␯. 共13兲 This decoupling is similar to the decoupling of the two- electron density matrix in the molecular spin-orbital basis. It may therefore be argued that⌬ and⌫ should be called the AO one- and two-electron density matrices, respectively.

However, in keeping with standard nomenclature, we shall refer toD as the AO density matrix, and to⌬and⌫as the matrices of expectation values of the one- and two-electron operators, respectively.

From the symmetry, trace, and idempotency properties of the one-electron density matrix,^26,27

D^†=D, 共14兲

TrDS=N_el, 共15兲

DSD=D, 共16兲

we straightforwardly obtain the following relations for⌬:

⌬^†=⌬, 共17兲

Tr⌬S⁻¹=N_el, 共18兲

⌬S⁻¹⌬=⌬. 共19兲 Although Eqs.共15兲and共16兲are formally equivalent to Eqs.

共18兲and共19兲, the relations for the standard AO density ma- trixDare somewhat simpler to use as they contain the metric S as opposed to the relations for ⌬, which contain the in- verted metricS⁻¹. We note that Eqs.共14兲–共19兲are necessary and sufficient conditions for a density matrix to represent a normalized single-determinant wave function.

B. Transformations of the density matrix

Next, we consider how ⌬ changes for a corresponding change in the reference wave function 兩0典 in Eq. 共5兲 de- scribed by the exponential operator

Tˆ= exp共i␬^ˆ兲, 共20兲

where␬^ˆ is a Hermitian one-electron operator,

␬^ˆ=

兺

␮␯ ␬␮␯a_␮^†a_␯, 共21兲

and␬is a Hermitian matrix. The transformed wave function may be expressed as

兩0˜典= exp共i␬^ˆ兲兩0典= exp共i␬^ˆ兲ai†a_j^†¯a_l^†兩vac典

=a_i^†a^†_j¯˜a_l^†兩vac典, 共22兲

where we have introduced the transformed creation operators

˜a_␮^† = exp共i␬^ˆ兲a_␮^† exp共−i␬^ˆ兲 共23兲 which satisfy the same anticommutation relations as the un- transformed operators

关˜a_␮^†,a˜_␯兴+=关exp共i␬^ˆ兲a_␮^†exp共−i␬^ˆ兲,exp共i␬^ˆ兲a_␯exp共−i␬^ˆ兲兴+

= exp共i␬^ˆ兲关a_␮^†,a_␯兴+exp共−i␬^ˆ兲=S_␯␮. 共24兲 The exponential operators of Eq.共20兲therefore represent the manifold of operators that conserve the metricS.

Using the Hausdorff expansion²⁶ and the anticommuta- tion relation of Eq. 共1兲, the transformed creation operators can be related to the untransformed ones as

a

˜_␮^† =a_␮^†+i关␬^ˆ,a_␮^†兴−1

2关␬^ˆ,关␬^ˆ,a_␮^†兴兴+ ¯

=a_␮^†+i

兺

␯ 共␬S兲_␮␯a_␯^†−1 2

兺

␯ 共␬S兲_␮␯² a_␯^†+ ¯

=

兺

␯ exp共i␬S兲_␮␯a_␯^†. 共25兲

In the special case whereS=I, exp共i␬兲 represents a unitary transformation of the orbitals.

The expectation value ofa_␮^†a_␯for the transformed state becomes

⌬

˜_␮␯=具0˜兩an†a_␯兩0˜典

=具0兩exp共−i␬^ˆ兲a_␮^†exp共i␬^ˆ兲exp共−i␬^ˆ兲a_␯exp共i␬^ˆ兲兩0典.

共26兲 Using Eqs.共23兲and共25兲, we straightforwardly obtain

exp共−i␬^ˆ兲a_␮^† exp共i␬^ˆ兲=

兺

exp共−i␬^ˆ兲a_␯exp共i␬^ˆ兲=

兺

␴ exp共iS␬兲_␯␴a_␴. 共27b兲 The insertion of these expressions into Eq.共26兲gives

⌬

˜ = exp共−iS^T␬^T兲⌬exp共i␬^TST兲, 共28兲 which describes how ⌬ transforms for a transformation of the reference state. If⌬fulfills Eqs.共17兲–共19兲, then so does

⌬˜ defined in Eq. 共26兲, as demonstrated in Appendix B. We conclude that⌬˜ fulfills Eqs.共17兲–共19兲and that exp共i␬^ˆ兲兩0典is a normalized single-determinant wave function. It may also be shown that all matrices fulfilling Eqs. 共17兲–共19兲 may be obtained by an appropriate choice of␬, so the transformation of Eq.共22兲constitutes a complete parametrization. Note that the matrix ␬ should comply with the projection relation ␬

=P共␬兲=P␬Q^†+Q␬P^†, whereP=DSandQ=I−DS, so as to eliminate redundancies—that is, sets of nonvanishing param- eters for which⌬˜共␬兲=⌬.

III. TIME-DEPENDENT EQUATIONS IN SECOND QUANTIZATION

A. Time-dependent wave functions and Hamiltonian The electronic system is described by the Hamiltonian Hˆ =Hˆ

0+Vˆ共t兲+Wˆ共t兲, 共29兲

where Hˆ

0 is the unperturbed operators, whereas Vˆ共t兲 and Wˆ共t兲 are first- and second-order perturbations that describe how the system interacts with the applied field. Att= −⬁, no field is present and no perturbations appear in the Hamil- tonian. Fort⬎−⬁, the fields are slowly applied according to the adiabatic approximation

Vˆ共t兲=

冕

Wˆ共t兲=

冕

冕

⫻exp兵关−i共␻1+␻2兲+ 2⑀兴t其d␻1d␻2, 共30b兲 where ⑀ is a real positive infinitesimal that ensures Vˆ共t= −⬁兲= 0. The perturbations are required to be Hermitian, so we have the relations

V˜共␻兲^†=V˜共−␻兲, 共31a兲

Wˆ共␻1,␻2兲^†=Wˆ共−␻1,−␻2兲. 共31b兲 We further require that Wˆ is symmetric in the two frequen- cies:

Wˆ共␻1,␻2兲=Wˆ共␻2,␻1兲. 共32兲 In response to these perturbations, the wave function changes,

兩0˜典= exp共i␬^ˆ共t兲兲兩0典, 共33兲

and we shall use this dependence to calculate the time de- pendence of expectation values induced by the perturbations.

We assume that the unperturbed system is described by an optimized single-determinant wave function 兩0典, for which the first-order variation of the energy具0兩H0兩0典 vanishes. Ex- pansion of the expectation value of the Hamiltonian operator as

具0˜兩Hˆ

0兩0˜典=具0兩Hˆ

0兩0典−i具0兩关␬^ˆ,Hˆ0兴兩0典+i²

2具0兩关␬^ˆ,关␬^ˆ,Hˆ0兴兴兩0典

− i³

3!具0兩关␬^ˆ,关␬^ˆ,关␬^ˆ,Hˆ0兴兴兴兩0典+O共␬^ˆ4兲 共34兲 then yields that the optimized state satisfies the stationary condition

d d␬␮␯具0˜兩H˜

0兩0˜典=i具0兩关Hˆ

0,a_␮^†a_␯兴兩0典= 0. 共35兲 We shall now consider how expectation values are affected by time-dependent perturbations.

B. The linear and quadratic response functions The linear and quadratic response functions describe the first- and second-order correction to the expectation value of a given 共time-independent兲 operator Aˆ when the system is subjected to a perturbation. To determine the linear and qua- dratic response functions, we begin by considering the time dependence of the expectation value 具0˜兩Aˆ兩0˜典 of a one- electron operator Aˆ, noting that, for this purpose, we only need to expand the wave function 兩0˜典of Eq. 共22兲to second order in the external perturbation:

␬^ˆ共t兲=␬^ˆ共1兲共t兲+␬^ˆ共2兲共t兲+ ¯ . 共36兲 The zero-order term␬^ˆ共0兲does not contribute since the unper- turbed wave function 兩0典satisfies the stationarity conditions Eq. 共35兲. To second order, we obtain the expansion

具0˜兩Aˆ兩0˜典=具0兩Aˆ兩0典−i具0兩关␬^ˆ共1兲共t兲,Aˆ兴兩0典

−¹₂具0兩关␬^ˆ共1兲共t兲,关␬^ˆ共1兲共t兲,Aˆ兴兴兩0典

−i具0兩关␬^ˆ共2兲共t兲,Aˆ兴兩0典. 共37兲 Since the response functions are defined in the frequency rather than in the time domain, we introduce the wave- function corrections in the frequency domain. By analogy with Eq.共30兲, we write

␬^ˆ共1兲共t兲=

冕

␬^ˆ共2兲共t兲=

冕冕

⫻exp兵关−i共␻1+␻2兲+ 2⑀兴t其d␻1d␻2, 共38b兲 where we require the second-order corrections to be symmet- ric in the frequencies:

␬^ˆ共2兲共␻1,␻2兲=␬^ˆ共2兲共␻2,␻1兲. 共39兲 Inserting the frequency expansions of the wave-function cor- rections of Eq.共38兲into Eq.共37兲, we obtain

具0˜兩Aˆ兩0˜典=具0兩Aˆ兩0典−i

冕

⬁

具0兩关␬^ˆ共1兲共␻兲,Aˆ兴兩0典

⫻exp关共−i␻⁺⑀兲t兴d␻

−1

2

冕冕

⫻exp关共−i共␻1+␻2兲+ 2⑀兲t兴d␻1d␻2

−i

冕冕

⫻exp兵关−i共␻1+␻2兲+ 2⑀兴t其d␻1d␻2. 共40兲 The integrand in the double commutator expression of Eq.

共40兲may be symmetrized in the integration variables␻1and

␻2using the operatorP₁₂, which creates the different permu- tations of the frequencies␻1 and␻2:

P₁₂具0兩关␬^ˆ共1兲共␻1兲,关␬^ˆ共1兲共␻2兲,Aˆ兴兴兩0典

=¹₂具0兩关␬^ˆ共1兲共␻1兲,关␬^ˆ共1兲共␻2兲,Aˆ兴兴兩0典 +¹₂具0兩关␬^ˆ共1兲共␻2兲,关␬^ˆ共1兲共␻1兲,Aˆ兴兴兩0典.

Comparing Eq.共40兲with the formal expansion of an expec- tation value in terms of response functions,¹

具0˜兩Aˆ兩0˜典=具0兩Aˆ兩0典+

冕

+

冕冕

⫻exp兵关−i共␻1+␻2兲+ 2⑀兴t其d␻1d␻2

+1 2

冕冕

⬁

具具Aˆ;Vˆ共␻1兲,Vˆ共␻2兲典典_␻1,␻2

⫻exp兵关−i共␻1+␻2兲+ 2⑀兴t其d␻1d␻2, 共41兲 we may identify the linear response function,

具具Aˆ;Vˆ共␻兲典典_␻= −i具0兩关␬^ˆ共1兲共␻兲,Aˆ兴兩0典, 共42兲 and the sum of a linear response function and a quadratic response function,

具具Aˆ;Wˆ共␻1,␻2兲典典_␻1+␻2+¹₂具具Aˆ;Vˆ共␻1兲,Vˆ共␻2兲典典_␻1,␻2

= −¹₂P₁₂具0兩关␬^ˆ共1兲共␻1兲,关␬^ˆ共1兲共␻2兲,Aˆ兴兴兩0典

−i具0兩关␬^ˆ共2兲共␻1,␻2兲,Aˆ兴兩0典. 共43兲 We shall later see how the linear and quadratic response functions in Eq.共43兲may be separated. The response of the system to a given perturbation can be determined from the wave-function correction of␬^ˆ. In the next section, the first- and second-order corrections are derived.

IV. FORM AND SOLUTION OF THE TIME-DEPENDENT EQUATIONS

A. Equations for the time development of the reference state

In this section, we first discuss the equations for deter- mining the time-dependent parameters ␬共t兲 and then derive the equations for determining the first- and second-order cor- rections of Eq.共36兲. We begin by rewriting␬^ˆ in component form,

␬^ˆ=_␮⬎␯

兺

兺

where the single-excitation operatorsa_␮^†a_␯have been divided into a set of AO excitations␮⬎␯ and a set of AO deexcita- tions␮⬍␯. As the AO excitations and deexcitations are for- mally equivalent, this division has no physical significance but will prove important when the paired structure of the response equations is investigated. Note that, in the AO rep- resentation, it is not possible to omit the diagonal operators a_␮^†a_␮共the number operators兲as done in the molecular-orbital 共MO兲 representation, where these operators may be viewed

as providing phase factors to be absorbed in the MO coefficients.²In Ref.11, the diagonal operators were mistak- enly assumed to vanish although, when the code was written,¹² no such assumption was made and the working equations and results are therefore valid.

For the further development, it is convenient to collect the operators of ␬^ˆ in a vector:

⌳=

冢

冣

where the three classes of operators are defined as

Q_I^†=a_␮^†a_␯, ␮⬎␯^, 共46兲

DI=a_␮^†a_␮, 共47兲

Q_I=a_␯^†a_␮, ␮⬎␯^. 共48兲 The parameters of ␬may similarly be arranged in a vector

␣^共i兲=

冢

冣

such that ␬^ˆ共i兲=␣^共i兲†·⌳=⌳^†·␣^共i兲. The superscript 共i兲 here in- dicates the order of the external perturbation, see Eq. 共36兲. The diagonal elements␬_␮␮have only a real component due to the Hermiticity of ␬^.

Likewise, to study the time evolution of兩˜0典 in the pres- ence of the time-dependent perturbation, we introduce the transformed operator basis

⌳˜ =

冢

冣

where Q˜

I= exp共i␬共t兲兲QIexp共−i␬共t兲兲 共51兲 and similarly forQ˜

I

†,D˜

I

†. The time evolution of兩0˜典may now be determined by invoking Ehrenfest’s theorem for the trans- formed operators of ⌳˜ ,

id

dt具0˜兩⌳˜兩0˜典=i具0˜兩

冉

冊

共52兲 This choice of operator simplifies the Ehrenfest theorem since the expectation value of⌳˜ is time independent:

d

dt具0˜兩⌳˜兩0˜典= d

dt具0兩⌳兩0典= 0 共53兲

yielding

−i具0˜兩

冉

冊

The Ehrenfest theorem results in a set of nonlinear equations that determine the time evolution of␬^˜共t兲.

In the next section, we solve the Ehrenfest equations order by order in the parameters ␬6共t兲. However, we first make a comment on notation. The supermatrix notation in- troduced above for collecting all AO excitation operators into a vector is essential for simplifying our subsequent deri- vations. In the remainder of this paper, we adopt Italic bold- face type to represent vectors and matrices in the supermatrix notation, keeping Roman boldface type to describe vectors and matrices in conventional matrix-vector notation. For ex- ample,b andcrepresent two vectors in the supermatrix no- tation, whereasbandcare conventional matrices of dimen- sion N², where N is the number of AOs: b^Tc=兺IbIcI

=兺_␮␯b_␮␯c_␮␯= Tr共b^Tc兲.

B. Expansion of the time-dependent equations in the external perturbation

To obtain an order-by-order expansion of Eq.共54兲, we first recognize that the time-derivative term in Eq.共54兲may be expressed as

i具0˜兩

冉

冊

+O共␬^ˆ3兲 共55兲

as shown in Appendix C. Expanding Eq.共54兲in orders of the external perturbation, using Eq. 共36兲, and collecting the terms linear in the perturbation, we obtain the first-order time-dependent equations

i具0兩关⌳,␬^ˆ˙共1兲兴兩0典−具0兩关⌳,关Hˆ

0,␬^ˆ共1兲兴兴兩0典= −i具0兩关⌳,Vˆ共t兲兴兩0典. 共56兲 Likewise, by collecting the second-order terms, we obtain the second-order time-dependent equations

i

2具0兩关⌳,关␬^ˆ共1兲,␬^ˆ˙共1兲兴兴兩0典−具0兩关⌳,␬^ˆ˙共2兲兴兩0典

=i具0兩关⌳,关Hˆ

0,␬^ˆ共2兲兴兴兩0典−1

2具0兩关⌳,关关Hˆ

0,␬^ˆ共1兲兴,␬^ˆ共1兲兴兴兩0典 +具0兩关⌳,Wˆ共t兲兴兩0典+i具0兩关⌳,关Vˆ共t兲,␬^ˆ共1兲兴兴兩0典. 共57兲

C. The first-order equations

The time-dependent equations, Eqs.共56兲and 共57兲, may be straightforwardly solved in frequency space. Using Eqs.

共30兲and共38兲, we obtain the first-order equation

冕

−具0兩关⌳,关Hˆ

0,␬^ˆ共1兲共␻兲兴兴兩0典兲d␻

= −

冕

⬁

exp关共−i␻⁺⑀兲t兴共i具0兩关⌳,Vˆ共␻兲兴兩0典兲d␻^, 共58兲 which gives the first-order response equations in the fre- quency space,

␻具0兩关⌳,␬^ˆ共1兲共␻兲兴兩0典−具0兩关⌳,关H0,␬^ˆ共1兲共␻兲兴兴兩0典

= −i具0兩关⌳,Vˆ共␻兲兴兩0典. 共59兲 Inserting␬^共1兲=⌳^†·␣^共1兲, we obtain

共具0兩关⌳,关Hˆ

0,⌳^†兴兴兩0典−␻具0兩关⌳,⌳^†兴兩0典兲·␣^共1兲共␻兲

=i具0兩关⌳,Vˆ共␻兲兴兩0典. 共60兲 Introducing the generalized Hessian and metric matrices

E_IJ^关2兴=具0兩关⌳I,关Hˆ

0,⌳J

†兴兴兩0典, 共61兲

S_IJ^关2兴=具0兩关⌳I,⌳J

†兴兩0典, 共62兲

and the property gradient vector

V_I^关1兴共␻兲=具0兩关⌳I,Vˆ共␻兲兴兩0典, 共63兲

the first-order response equations关Eq. 共60兲兴may, in 共super兲 matrix notation, be expressed as

共E^关2兴−␻^S关2兴兲␣^{共1兲=iV关1兴}共␻兲. 共64兲

Note that E^关2兴 is the second derivative of the energy with respect to the electronic-structure parameters 关see the second-order term in Eq. 共34兲兴and thus positive definite for an optimized ground state.

D. The second-order equations

Invoking Eqs.共30兲and共38兲, we may express the second- order equations in the frequency domain as

冕冕

⬁

e^{共−i␻1+⑀兲t}e^{共−i␻2+⑀兲t}

冉

冊

=

冕冕

冉

+i具0兩关⌳,关Vˆ共␻2兲,␬^ˆ共1兲共␻1兲兴兴兩0典−i共␻1+␻2兲具0兩关⌳,␬^ˆ共2兲共␻1,␻2兲兴兩0典

冊

Using the definitions ofE^关2兴andS^关2兴from Eqs.共61兲and共62兲, together with␬^ˆ共2兲共␻1,␻2兲=⌳^†·␣^共2兲共␻1,␻2兲, we obtain

兵E^关2兴−共␻1+␻2兲S^关2兴其␣^共2兲共␻1,␻2兲

=P₁₂

冋

+ i

2␻1具0兩关⌳,关␬^ˆ共1兲共␻1兲,␬^ˆ共1兲共␻2兲兴兴兩0典

−具0兩关⌳,关Vˆ共␻1兲,␬^ˆ共1兲共␻2兲兴兴兩0典

+i具0兩关⌳,Wˆ共␻1,␻2兲兴兩0典

册

In Eq.共66兲, we have performed a symmetrization of the in- tegration variables to obtain response parameters ␬^共2兲

⫻共␻1,␻2兲 that are symmetric in the frequency indices.

The second-order equations may be expressed in a more compact form, inserting␬^ˆ共1兲共␻兲=⌳^†·␣^共1兲共␻兲and introducing the three-index supermatricesE^关3兴 andS^关3兴,

E_IJK^关3兴 =¹₂具0兩关⌳I,关关Hˆ

0,⌳J†兴,⌳K†兴兴兩0典, 共67兲

S_IJK^关3兴 =¹₂具0兩关⌳I,关⌳J†,⌳K†兴兴兩0典. 共68兲 We furthermore introduceV^关2兴共␻兲, obtained fromE^关2兴 in Eq.

共61兲 by replacingHˆ

0 withVˆ共␻兲, and W^关1兴共␻1,␻2兲 obtained by replacing Vˆ共␻兲 by Wˆ共␻1,␻2兲 in Eq. 共63兲. The second- order equations can then be written as

兺

=P₁₂

兺

⫻共␻1兲␣J共1兲共␻2兲+iW_I^关1兴共␻1,␻2兲兴. 共69兲

V. RESPONSE FUNCTIONS AND THEIR RESIDUES A. The linear response functions

The computational expression for the linear response function may be obtained by inserting the first-order wave- function correction关Eq. 共64兲兴into the linear response func- tion in Eq.共42兲. Renaming the perturbation operatorVˆ共␻兲as Bˆ and introducing

A_J^关1兴= −具0兩关⌳J†,Aˆ兴兩0典, 共70a兲

B_J^关1兴=具0兩关⌳J,Bˆ兴兩0典, 共70b兲

we obtain

具具Aˆ;Bˆ典典_␻= −A^关1兴共E^关2兴−␻^S关2兴兲⁻¹B^关1兴. 共71兲 The linear response function may thus be calculated by solv- ing one set of linear equations at each frequency,

N^b共␻兲=共E^关2兴−␻^S关2兴兲⁻¹B^关1兴. 共72兲 B. Pairing of the excitation energies obtained from the linear response function

The excitation energies are identified as the poles of the linear response function of Eq.共71兲. The excitation energies are therefore solutions to the generalized eigenvalue problem E^关2兴Xf=S^关2兴Xf␻f. 共73兲 In the standard formulation of response theory in the MO basis, it is has been shown that the excitation energies are paired—that is, if ␻f is an eigenvalue, then −␻f is also an eigenvalue. The standard MO proof of pairing² cannot be applied directly in the AO basis due to the presence of the diagonal operators D_I. We first analyze here the structure of E^关2兴andS^关2兴and then show how the pairing is obtained in the AO basis. Note that when pairing was discussed in Ref.11, consideration was not given to the diagonal operators DI.

Partitioning of⌳into the three classes 关Eq.共46兲–共48兲兴, we may writeE^关2兴as

E^关2兴=

冢

0,D兴兴兩0典 具0兩关Q,关Hˆ

0,Q兴兴兩0典 具0兩关D,关Hˆ

0,Q^†兴兴兩0典 具0兩关D,关Hˆ

0,D兴兴兩0典 具0兩关D,关Hˆ

0,Q兴兴兩0典 具0兩关Q^†,关Hˆ

0,Q^†兴兴兩0典 具0兩关Q^†,关Hˆ

0,D兴兴兩0典 具0兩关Q^†,关Hˆ

0,Q兴兴兩0典

冣

The nine blocks in Eq.共74兲can all be written in terms of the following four matrices:

AIJ=具0兩关QI,关Hˆ

0,Q_J^†兴兴兩0典, 共75兲

BIJ=具0兩关QI,关Hˆ

0,QJ兴兴兩0典, 共76兲

LIJ=具0兩关QI,关Hˆ

0,DJ兴兴兩0典, 共77兲

NIJ=具0兩关QI,关Hˆ

0,DJ兴兴兩0典. 共78兲

For example, the block具0兩关Q^†,关Hˆ

0,Q兴兴兩0典can be written in terms ofA, exploiting the fact that any scalar can be written as A=共A^†兲^*,

具0兩关QI

†,关Hˆ

0,QJ兴兴兩0典=共具0兩关QI

†,关Hˆ

0,QJ兴兴兩0典^†兲^*

=具0兩关Q_I,关Hˆ

0,Q_J^†兴兴兩0典^*=A_IJ^*. 共79兲 To rewrite the blocks 具0兩关D,关Hˆ

0,Q^†兴兴兩0典 and 具0兩关D,关Hˆ

0,Q兴兴兩0典, the Jacobi identity关see Eq.共C9兲兴and the stationarity condition Eq.共35兲are used,

具0兩关DI,关Hˆ

0,Q_J^†兴兴兩0典

= −具0兩关Hˆ

0,关QJ

†,DI兴兴兩0典−具0兩关QJ

†,关DI,Hˆ

0兴兴兩0典

=具0兩关Q_J^†,关Hˆ

0,D_I兴兴兩0典

=共具0兩关QJ

†,关Hˆ

0,DI兴兴兩0典^†兲^*

=共具0兩关QJ,关Hˆ

0,DI兴兴兩0典兲^*=L_JI^*, 共80兲

giving 具0兩关D,关Hˆ

0,Q^†兴兴兩0典=L^†. 共81兲

Rewriting the remaining blocks in a similar way, we obtain

E^关2兴=

冢

冣

The matrixS^关2兴 may, in a similar way, be written as

S^关2兴=

冢

冣

where

⌺IJ=具0兩关QI,Q_J^†兴兩0典, 共84兲

⌬IJ=具0兩关Q_I,Q_J兴兩0典, 共85兲

TIJ=具0兩关QI,DJ兴兩0典, 共86兲

⑀IJ=具0兩关DI,DJ兴兩0典. 共87兲

The block⑀ containing two diagonal operators vanishes for real orbitals,

⑀IJ=具0兩关DI,DJ兴兩0典具0兩关a_␮^†a_␮,a_␯^†a_␯兴兩0典

=S_␮␯具0兩a_␮^†a_␯兩0典−S_␯␮具0兩a_␯^†a_␮兩0典= 0. 共88兲 Note thatE^关2兴 andS^关2兴 are both Hermitian. ForE^关2兴, Hermi- ticity follows from the fact thatA^†=A,B^T=B, andNis real and symmetric. For S^关2兴, it occurs since ⌺^†=⌺, ⌬^†= −⌬^*, and⑀is imaginary and antisymmetric.

Let us now assume that the vector

X_f=

冢

冣

is an eigenvector for Eq.共73兲with eigenvalue␻f,

冢

冣 冢

冣

冢

冣 冢

冣

共90兲 We may now express Eq. 共90兲 in component form as three sets of equations,

AZ+LU+BY=␻f共⌺Z+␪U+⌬Y兲, 共91a兲

L^†Z+NU+L^TY=␻f共␪^†Z+⑀^U−␪^TY兲, 共91b兲 B^*Z+L^*U+A^*Y=␻f共−⌬^*Z−␪^*U−⌺^*Y兲. 共91c兲 We will now prove that the paired vectorX_−f given as

X_−f=

冢

冣

is an eigenvector for Eq.共73兲with eigenvalue −␻f,

冢

冣 冢

冣

= −␻f

冢

冣 冢

冣

In component form, Eq.共93兲becomes

AY^*+LU^*+BZ^*= −␻f共⌺Y^*+␪^U*+⌬Z^*兲, 共94a兲 L^†Y^*+NU^*+L^TZ^*= −␻f共␪^†Y*−⑀^U*−␪^TZ*兲, 共94b兲 B^*Y^*+L^*U^*+A^*Z^*= −␻f共−⌬^*Y^*−␪^*U*−⌺^*Z^*兲.

共94c兲 Complex conjugation of Eq. 共94兲reveals the equivalence of Eqs. 共91a兲 and共94c兲, of Eqs. 共91b兲 and共94b兲, and of Eqs.

共91c兲 and共94a兲. Therefore, if X_f is an eigenvector with ei- genvalue␻f,

E^关2兴Xf=␻fS^关2兴Xf, 共95兲 thenX_−f is an eigenvector with eigenvalue −␻f,

E^关2兴X_−f= −␻fS^关2兴X_−f. 共96兲

C. The residue of the linear response functions In the previous section, the solution vectors to the gen- eralized eigenvalue problem关Eq.共73兲兴were shown to appear in pairs and Eqs. 共96兲and共95兲may be combined as

1

␻f

E^关2兴X_f= sgn共f兲S^关2兴X_f. 共97兲