TDHF + CIS in Python

So you completed a Hartree-Fock procedure, and you even transformed your two electron integrals. Now what can you do? We can use those results, specifically the orbital energies (eigenvalues of the Fock matrix) and the transformed two-electron integrals (using the eigenvalues of the Fock matrix), to calculate some simple response properties, such as excitation energies!. Configuration Interaction Singles (CIS), is a subset of the TDHF equations. If you haven’t performed the previous calculations, I’ll include my initial eigenvalues and transformed two electron integrals at the end. The values are for HeH+ at a bond length of 0.9295 Angstrom with an STO-3G basis set. We create and solve:

With the elements of \({\mathbf{A}}\) and \({\mathbf{B}}\) as

and

Now, when we transformed our two electron integrals from atomic orbital (AO) basis to molecular orbital (MO) basis, we made the implicit assumption that we were working under a closed-shell approximation. Electrons have spin, and we assumed we had an even number of electrons so that all their spins paired and integrated out. If you’ve taken any general chemistry, you know that we fill each orbital with two electrons: one spin up, and one spin down. (Electrons are fermions, which means two electrons of opposite spin can occupy the same spatial location. Their opposite, bosons, are a little weirder, which is why we can get Bose-Einstein condensates — think a lot of particles occupying the same place in space). Anyway, we make the assumption that opposite spin electrons can share the same spatial orbital, which is reasonable for most calculations, but not strictly true – after all, the definition of an orbital is a one electron wavefunction. Moving along now…

We need to explicitly account for the spin of the electrons, which means we need to satisfy this transformation:

Where \(\omega\) is the spin of the electron, with coordinates 1 and 2. Note that \((pq\vert rs) \neq \langle pq \vert rs\rangle\), this is because we will be moving from chemists’ notation to physicists’ notation, which dominates the literature. Since the CIS and TDHF equations need only double bar integrals, e.g.

or, showing the conversion between physicists’ and chemists’ notation more directly

We can also account for that in our transformation. In Python this gives:

# Create spin-adapted two electron integral  
# in double-bar format in a 4D array 'ints'  
# teimo(i,j,k,l) grabs MO basis (pq | rs) value  
# Note index change from starting at '1' to '0'  
ints=np.zeros((dim*2,dim*2,dim*2,dim*2))  
for p in range(1,dim*2+1):  
    for q in range(1,dim*2+1):  
        for r in range(1,dim*2+1):  
            for s in range(1,dim*2+1):  
                value1 = teimo((p+1)//2,(r+1)//2,(q+1)//2,(s+1)//2) * (p%2 == r%2) * (q%2 == s%2)  
                value2 = teimo((p+1)//2,(s+1)//2,(q+1)//2,(r+1)//2) * (p%2 == s%2) * (q%2 == r%2)  
                ints[p-1,q-1,r-1,s-1] = value1 - value2

Now, because Python begins its indexing at 0, and I began my indexing at 1, I had to make the index change at the end. From now on, the indexing starts at 0, so that \((11\vert 11) = (00\vert 00)\) and so on. Now that we have our integrals, we can also spin adapt our orbital energies. This basically maps spatial orbital energy 1 (which contains two electrons) to spin orbital 1 and 2: odd numbers are spin up, even numbers are spin down. If the spatial orbitals are found in an array E, now moved to an array of double the dimension fs:

 
# Spin basis fock matrix eigenvalues  
fs = np.zeros((dim*2))  
for i in range(0,dim*2):  
    fs[i] = E[i//2]  
fs = np.diag(fs)  

Simple enough.

Putting it altogether then (I hard coded the initial values for self-containment):

 
#!/usr/bin/python

####################################
#
#  CI Singles (CIS) and
#
#  Time-Dependent Hartree Fock (TDHF)
#
####################################

from __future__ import division
from __future__ import print_function
import math
import numpy as np

####################################
#
#   FUNCTIONS
#
####################################

# Return compound index given four indices
def eint(a,b,c,d):
  if a > b: ab = a*(a+1)/2 + b
  else: ab = b*(b+1)/2 + a
  if c > d: cd = c*(c+1)/2 + d
  else: cd = d*(d+1)/2 + c
  if ab > cd: abcd = ab*(ab+1)/2 + cd
  else: abcd = cd*(cd+1)/2 + ab
  return abcd

# Return Value of spatial MO two electron integral
# Example: (12|34) = tei(1,2,3,4)
def teimo(a,b,c,d):
  return ttmo.get(eint(a,b,c,d),0.0e0)

####################################
#
#  INITIALIZE ORBITAL ENERGIES 
#  AND TRANSFORMED TWO ELECTRON
#  INTEGRALS  
#
####################################

Nelec = 2 # we have 2 electrons in HeH+
dim = 2 # we have two spatial basis functions in STO-3G
E = [-1.52378656, -0.26763148]
ttmo = {5.0: 0.94542695583037617, 12.0: 0.17535895381500544, 14.0: 0.12682234020148653, 17.0: 0.59855327701641903, 19.0: -0.056821143621433257, 20.0: 0.74715464784363106}

####################################################
#
#  CONVERT SPATIAL TO SPIN ORBITAL MO
#
####################################################

# This makes the spin basis double bar integral (physicists' notation)

spinints=np.zeros((dim*2,dim*2,dim*2,dim*2))
for p in range(1,dim*2+1):
  for q in range(1,dim*2+1):
    for r in range(1,dim*2+1):
      for s in range(1,dim*2+1):
        value1 = teimo((p+1)//2,(r+1)//2,(q+1)//2,(s+1)//2) * (p%2 == r%2) * (q%2 == s%2)
        value2 = teimo((p+1)//2,(s+1)//2,(q+1)//2,(r+1)//2) * (p%2 == s%2) * (q%2 == r%2)
        spinints[p-1,q-1,r-1,s-1] = value1 - value2

spatints=np.zeros((dim,dim,dim,dim))
for p in range(1,dim+1):
  for q in range(1,dim+1):
    for r in range(1,dim+1):
      for s in range(1,dim+1):
        spatints[p-1,q-1,r-1,s-1] = teimo(p,r,q,s)

#####################################################
#
#  Spin basis fock matrix eigenvalues 
#
#####################################################

fs = np.zeros((dim*2))
for i in range(0,dim*2):
    fs[i] = E[i//2]
fs = np.diag(fs)

######################################################
#
#   CIS & TDHF CALCULATION 
#
######################################################

NOV = Nelec*(2*dim-Nelec)

A = np.zeros((NOV,NOV))
B = np.zeros((NOV,NOV))
I = -1
for i in range(0,Nelec):
  for a in range(Nelec,dim*2):
    I = I + 1
    J = -1
    for j in range(0,Nelec):
      for b in range(Nelec,dim*2):
        J = J+1
        A[I,J] = (fs[a,a] - fs[i,i]) * ( i == j ) * (a == b) + spinints[a,j,i,b]
        B[I,J] =  spinints[a,b,i,j]


# Solve CIS matrix equation
ECIS,CCIS = np.linalg.eig(A)
print("E(CIS) = ", np.amax(ECIS))
# Solve TDHF matrix equation
M = np.bmat([[A,B],[-B,-A]])
ETD,CTD = np.linalg.eig(M)
print("E(TDHF) = ", abs(np.amax(ETD)))

Running the code we find our first excitation energy at the CIS level is E(CIS) = 0.911, and our first excitation at the TDHF level is E(TDHF) = 0.902, (all in Hartrees) in agreement with the literature values in the paper here. A couple of notes about the method. First, the TDHF method can be seen as an extension of the CIS approach. The \(A\) matrix in the TDHF working equations is just the CIS matrix. This also means that the TDHF matrix is twice the dimension of the CIS matrix ( four times as large), which can get very costly. With a little bit of (linear) algebra, we can get the TDHF equations in the form:

Which is the same dimension as the CIS matrix. This is why you rarely see CIS in practice: for the same cost, you can get a TDHF calculation, which is more accurate because it includes the block off-diagonal elements. I also want to mention that TDDFT takes almost exactly the same form. The only difference is that the replace the two electron integrals with the exchange elements in DFT, and the HF orbital energies are the Kohn-Sham(KS) orbital eigenvalues. Because DFT generally handles correlation better, most single excitations are modeled using TDDFT. Finally, we did solve for all the excitation energies in the minimal basis…but in general we don’t need them all, and it is way to costly to do so. In this case, we can use iterative solvers — like the Davidson iteration — to pick off just the lowest few eigenvalues. Gaussian09, by default, solves for just the lowest three eigenvalues unless you ask for more.

For your reference, the transformed two electron integrals \((pq\vert rs)\) are given below, first four columns are index of p,q,r,s, and the last column is the value:

1 1 1 1 0.94542695583  
1 1 1 2 0.175358953815  
1 1 2 1 0.175358953815  
1 1 2 2 0.598553277016  
1 2 1 1 0.175358953815  
1 2 1 2 0.126822340201  
1 2 2 1 0.126822340201  
1 2 2 2 -0.0568211436214  
2 1 1 1 0.175358953815  
2 1 1 2 0.126822340201  
2 1 2 1 0.126822340201  
2 1 2 2 -0.0568211436214  
2 2 1 1 0.598553277016  
2 2 1 2 -0.0568211436214  
2 2 2 1 -0.0568211436214  
2 2 2 2 0.747154647844

Addendum: Getting transition dipole moment

Let’s say you have your z-component dipole moment integrals (in the molecular orbital basis) collected into a matrix \(Z_{pq} = \langle p \vert \overrightarrow{z} \vert q \rangle\).

You can compute the z-component of the transition dipole moment from your transition density, which is either the appropriate eigenvector for a given excitation \(X_{ia}\) (CIS) or \(X_{ia}\) and \(Y_{ai}\) (TDHF). If the dipole integral matrix is \(N \times N\), you’ll want to reshape your \(X_{ia}\) and \(Y_{ai}\) to be \(N \times N\) as well, instead of \(OV \times 1\). This gives you a \(N \times N\) transition density matrix of

The expectation value of the z-component of this transition dipole moment is then the trace of the density dotted into the z-component dipole integrals, e.g.

Then do the same for \(\langle x \rangle\) and \(\langle y \rangle\) and sum each component for the total transition dipole.

You’ll also need to repeat the process for each transition you are interested in.

Note that most dipole moment integrals are computed in the atomic orbital basis, so don’t forget to transform your transition density to the AO basis (or dipole integrals to MO–it doesn’t matter which way, expectation values are independent of basis–the key is that both matrices are in the same basis).

Efficient Two-Electron Integral Transformations in Python, or, Adventures in Scaling

Low-scaling algorithms are important in any computational development, and this is never more true than in quantum chemistry. A golden example is in the integral transformation routines in electronic structure packages. Without a smarter algorithm, the transformations scale as \(O(N^8)\) (yikes!!), making computations on large systems nearly impossible. But with a little thought, we can reign in that transformation to a (not nearly as bad) \(N^5\) method.

If you have performed a Hartree-Fock (HF) calculation, you are left with a matrix of coefficients. Mathematically, these coefficients are the eigenvectors that correspond to your HF Hamiltonian (Fock matrix). The eigenvalues are your orbital energies. So far so good. The real accurate (and interesting) work in electronic structure theory comes afterward, where we either refine the energies we have obtained or calculate molecular properties like UV-Vis spectra and so on. Most of these methods – actually, all these methods – require transforming your two electron integrals from an atomic orbital basis (your standard run-of-the-mill basis functions are your AO basis) into a molecular orbital basis. Put another way, now that you have solved the wavefunction of the Hartree-Fock system, you need to project your two electron integrals onto that wavefunction. The way you do this looks like so:

In code, that looks something like:

for p in range(0,dim):  
    for q in range(0,dim):  
        for r in range(0,dim):  
            for s in range(0,dim):  
                for mu in range(0,dim):  
                    for nu in range(0,dim):  
                        for lam in range(0,dim):  
                            for sig in range(0,dim):  
                                TxInt[p,q,r,s] += C[p,mu]*C[q,nu]*C[r,lam]*C[s,sig]*UnTxInt[mu,nu,lam,sig]

Holy cow. EIGHT loops of dimension number of basis functions. That should scale on the order of \(N^8\). Few methods in electronic structure theory scale worse than that (though they do exist…here’s lookin’ at you Full CI). This means that for most calculations, if you use this method of integral transformation, this transformation will be your most expensive step. Thankfully, we can come up with a smarter way to do the above transformation. Because the coefficients are independent (i.e \(C^{p}_\mu\) is independent of \(C^{q}_\nu\)), we can rewrite our first equation as

Or, like this, which is a lot clearer to me:

Where we perform four “quarter-transformations”, save each transformation, and use in the next transformation. Doing it this way gives us four steps of 5-dimension loops, so it should scale on the order of \(N^5\). In code, that looks like:

for p in range(0,dim):  
    for mu in range(0,dim):  
        temp[p,:,:,:] += C[p,mu]*UnTXInt[mu,:,:,:]  
    for q in range(0,dim):  
        for nu in range(0,dim):  
            temp2[p,q,:,:] += C[q,nu]*temp[p,nu,:,:]  
        for r in range(0,dim):  
            for lam in range(0,dim):  
                temp3[p,q,r,:] += C[r,lam]*temp2[p,q,lam,:]  
            for s in range(0,dim):  
                for sig in range(0,dim):  
                    TxInt[p,q,r,s] += C[s,sig]*temp3[p,q,r,sig]

Much nicer. You’ll notice that we had to pre-allocate the ‘temp’ matrices, to store the results between quarter transformations. The transformation also makes use of the ‘slice’ notation in Python/ NumPy. Using this, we perform a transformation over a whole dimension, instead of one index at a time. It’s a little weird to be working with full dimensions, instead of just indices, but it works well. Here is the full code, with random integer arrays built in to act as our toy four-dimensional integrals. The toy matrix of coefficients, is, like all matrices, 2D. I built in a check, so you can compare the two methods. It spits out two transformed integrals with randomly chosen indices – if/when you run it, you should make sure that the values match. If they don’t, something is wrong!

 
#!/usr/bin/python

####################################  
#  
# SAMPLE INTEGRAL TRANSFORMATION CODE  
#  
####################################

from __future__ import division  
import sys  
import math  
import numpy as np  
import time

#####################################

# Initialize Arrays  
dim = 2 # dimension of arrays ... e.g number of basis functions  
MO1 = np.zeros((dim,dim,dim,dim)) # For our first dumb O[N^8] method  
MO2 = np.zeros((dim,dim,dim,dim)) # For our smarter O[N^5] method

INT = np.random.randint(9,size=(dim,dim,dim,dim)) # Our toy "two electron integrals"  
C = np.random.randint(9,size=(dim,dim)) # Toy "wavefunction coefficients"

# Begin first method. It scales as N^8, as you could  
# have guessed with there being 8 loops over dimension 'dim' (N)

t0 = time.time()  
for i in range(0,dim):  
    for j in range(0,dim):  
        for k in range(0,dim):  
            for l in range(0,dim):  
                for m in range(0,dim):  
                    for n in range(0,dim):  
                        for o in range(0,dim):  
                            for p in range(0,dim):  
                                MO1[i,j,k,l] += C[i,m]*C[j,n]*C[k,o]*C[l,p]*INT[m,n,o,p]  
t1 = time.time()

# Begin second method, scaling as N^5. We end up having four 5-loops, each  
# over dimension 'dim' (N).

t2 = time.time()  
temp = np.zeros((dim,dim,dim,dim))  
temp2 = np.zeros((dim,dim,dim,dim))  
temp3= np.zeros((dim,dim,dim,dim))  
for i in range(0,dim):  
    for m in range(0,dim):  
        temp[i,:,:,:] += C[i,m]*INT[m,:,:,:]  
    for j in range(0,dim):  
        for n in range(0,dim):  
            temp2[i,j,:,:] += C[j,n]*temp[i,n,:,:]  
        for k in range(0,dim):  
            for o in range(0,dim):  
                temp3[i,j,k,:] += C[k,o]*temp2[i,j,o,:]  
            for l in range(0,dim):  
                for p in range(0,dim):  
                    MO2[i,j,k,l] += C[l,p]*temp3[i,j,k,p]  
t3 = time.time()

# Set up random index to check correctness.  
i = np.random.randint(dim)  
j = np.random.randint(dim)  
k = np.random.randint(dim)  
l = np.random.randint(dim)

print MO1[i,j,k,l]  
print MO2[i,j,k,l]  
print "TIME1: ", t1-t0  
print "TIME2: ", t3-t2  

When I ran the code, moving from a dimension of 4 to a dimension of 8 (e.g I doubled the basis functions), the first method went from 0.29 seconds to 71.5 seconds, a jump of 246 times longer, versus the second method, which went from 0.01 seconds to 0.29 seconds, a jump of only 29 times. This is almost exactly as predicted. Doubling the basis for an \(N^8\) method gives \(2^8 = 256\) times longer, and doubling the basis for an \(N^5\) algorithm gives \(2^5 = 32\) times longer.

The new method is also very amenable to parallelization. Most electronic structure computations need to be parallelized, as model systems get larger and larger. Note that the \(N^5\) method is performed in four independent steps. Because of this independence, we can make our code run in parallel, and perform the quarter transformations on separate processors.

A great example why choosing your algorithm matters in quantum chemistry!

Awk Script for Gaussian TD-DFT/TD-HF Analysis

When I am generating UV-Vis spectra (in some cases near-IR spectra) for some system, I am often dealing with hundreds of individual electronic transitions. And most of the time, these transitions don’t even matter (e.g. their oscillator strength is ~ 0). This happens a ton when I am generating spectra for semiconductor quantum dots. TD analysis shows tons of ‘transitions’ but they exist in the band gap, and I’d rather filter them out. See, for example, below:

TD-DFT actually shows tons of transitions from ~2 to 4.5 eV…they just don’t matter! So for my own sanity in pouring through the data, I wrote this Awk script (well, Awk in a Bash script) to pull out the transitions I want, based on oscillator strength (intensity). If you don’t know Awk, you really should. It has nasty syntax (like Perl…), but once you get used to it, it makes text editing and parsing SO SO EASY. I love it. I use Awk all the time for my work. Take a look, and feel free to grab a copy for yourself. To run the script, just type (minus the quotes, and assuming you saved it as ‘tdanalysis’):

./tdanalysis [FILE] [minimum oscillator strength]

Enjoy!

#!/bin/bash

if ["$#" == "0"]; then  
 printf "********************\nTD_ANALYSIS\nPrints TDDFT Excitations above a given oscillator strength.\nUSE: td_analysis [LOG FILE] [MINIMUM OSCILLATOR STRENGTH]\n********************\n"  
 exit 1  
fi

# Grab file (Gaussian .log file) and desired minimum oscillator strength

# The MULT variable is used to give you

# % contributions based on whether the system is open or closed shell.

FILE=$1  
OSCIL=$2  
MULT=`gawk '/Charge/ && /Multiplicity/ {print $6%2+1}' $1`  
gawk -v mult=$MULT 'BEGIN {format = "%-6s %-10s %10s %10s %6s\n"  
 print "TD Analysis of",ARGV[1]}  
 /f=([0-9])/ { printf format, "\nEnergy [eV]:",$5, " ","Oscillator Strength:", substr($9,3,6)"\n ----------------------------------------------------------------"};  
 /->/ {printf format, "", $1$2$3,$4," ","("mult*100*($4)^2"%)"}  
 END { print "\nEnd of Analysis\n\n"}' $FILE \vert  gawk -v oscil=$OSCIL 'RS="\n\n" {if($6 > oscil) {print $0"\n"}}'

You may need to change gawk to awk, depending on your system. I’ve never had problems with gawk though.

Derivation of Time Dependent Hartree-Fock (TDHF) Equations

This derivation is largely based off the TD-DFT derivation (TDDFT and TDHF are very similar), found in the wonderful review, “Single-Reference ab Initio Methods for the Calculation of Excited States of Large Molecules” by Dreuw and Head-Gordon (Chem Reviews, 2005).

We start with the time-dependent Schrödinger equation, given in Hartree-Fock form, along with its adjoint:

Now,

Substituting the first two expressions into the right hand side, these two expressions yields the Dirac form of the time-dependent Hartree-Fock equations,

where \({\mathbf{F}}\) is the Fock matrix and \({\mathbf{P} = \mathbf{CC}^{\dagger}}\) is the density matrix, same as in the time independent Hartree-Fock equations. Before a time dependent perturbation is applied, we assume that the system is in its electronic ground state, such that

and

with \({\mathbf{P}^{(0)}}\) and \({\mathbf{F}^{(0)}}\) as the unperturbed density and Fock matrix, respectively. The first condition comes from the fact that commuting matrices share a common set of eigenvectors, and the second comes from the fact that at convergence the eigenvectors of the Fock matrix are orthonormal. (Recall that the density matrix may be written as \({\mathbf{P} = \mathbf{CC}^{\dagger}}\), where the columns of \({\mathbf{C}}\) form an orthonormal set, e.g. \({\mathbf{C}^{\dagger}\mathbf{C} = \mathbf{1}}\)). The Fock matrix elements are given by

Note the Fock matrix dependence on the density matrix. At convergence, the Fock matrix is diagonal, with the elements corresponding to the orbital energies \({\epsilon}\).

Furthermore, the density matrix enjoys the following relations:

Following convention, indices \({\{i,j,...\}}\) correspond to occupied orbitals, \({\{a,b,...\}}\) correspond to unoccupied orbitals, and \({\{p,q,...\}}\) correspond to general orbitals. In general perturbation theory, a perturbed wavefunction (or density matrix, for our purposes) can be decomposed, to first order, as

and likewise for the density-dependent Fock matrix

where the superscript (1) indicates the first-order time dependent change. Inserting the expressions 13 and 14 into the time-dependent equation 6 and collecting only the first order terms yields

The time-dependent perturbation can be described by a single Fourier component,

The matrix \({f}\) is a one-electron operator, which describes the applied perturbation. This perturbation acts on the electron density, resulting in a first order Fock matrix response

so that the overall first order change in the Fock matrix is

Similarly, the first order density response is given as

with \({d_{pq}}\) as the perturbation densities. Plugging the above expressions for \({\mathbf{F}^{(1)}}\) and \({\mathbf{P}^{(1)}}\) into the first order TDHF equation 15, and collecting the terms containing \({e^{-i\omega t}}\) gives (after some algebra):

The terms multiplied by \({e^{i \omega t}}\) give the complex conjugate of the above expression. Because the density matrix must be idempotent, we can show that

If \({\mathbf{PP} = \mathbf{P}}\), and we expand \({\mathbf{P}}\) as a perturbation with arbitrary scalar \({\lambda}\), we have:

Equating like powers of \({\lambda}\) yields the expression 21, as well as our original expression 8. Using 21 allows us to restrict the elements of \({d_{pq}}\), utilizing the nature of the zeroth order density matrix given in 11. As an example, consider the element \({d_{ii}}\):

which is only true if \({P_{ii}^{(1)} = 0}\). Similar arguments show that the only contributing blocks of \({d_{pq}}\) are \({d_{ia}}\) and \({d_{ai}}\). In other words, the only contributions to the TDHF equations are the occupied-virtual and virtual-occupied blocks. All virtual-virtual and occupied-occupied blocks are necessarily zero. With this in mind, along with the diagonal nature of the unperturbed Fock matrix, yields

and

with \({x_{ai} = d_{ai}}\) and \({y_{ai} = d_{ia}}\), as is convention. The derivatives inside the expression evaluate as follows:

Finally, we make the assumption that the perturbation is infinitely small, which sends \({f_{ai} = f_{ia} \rightarrow 0}\), and we recover (recognizing that \({F_{pp}^{(0)} = \epsilon_p}\) and \({P_{ii}^{(0)} = 1}\)):

Which is the non-Hermitian eigenvalue TDHF equation. The elements of \({\mathbf{A}}\) and \({\mathbf{B}}\) are

and

Alternate derivation of linear response function

I’ve given a derivation before of the density-density response function, but today I want to give an alternate derivation of the more general linear-response function, which will prove to be useful in the derivation of the time-dependent Hartree-Fock equations (TDHF), also known in the nuclear physics community as the Random Phase Approximation (RPA). This derivation is largely taken from McWeeny (1989).

In order to derive the response function — also called frequency-dependent polarizability — we must first partition the Hamiltonian into two parts: the time-independent Hamiltonian, and the time dependent response:

Furthermore, the time-dependent Schrodinger equation is given as:

In the interaction picture, \({\Psi'_0(t) = \Psi'_0e^{-iE_0t/\hbar}}\). Thus we can partition the time-dependent wavefunction (expanding over basis of complete eigenstates) as

Where

are real, positive, exact excitation frequencies of unperturbed system. Note that \({c_0(t) = 1}\). We are assuming that \({H'(t)}\) is turned on slowly at time \({t \rightarrow -\infty}\). Substitute 1 and 3 into 2, separate the orders, and impose the boundary conditions that \({c_0 = 1}\) and \({c_m = 0}\) \({(n\neq0)}\) at \({t \rightarrow -\infty}\). This gives

If we let

Where a ‘fixed’ Hermitian operator \({A}\) determines the `shape’ of the perturbation, while time dependence is confined to the (real) ‘strength’ factor \({F(t)}\). For a perturbation beginning at time \({t \rightarrow -\infty}\) up to time \({t}\),

Which, to first order, determines the perturbed wavefunction. Now we are interested not in the perturbed wavefunction per se, but rather in the response of an observable \({O}\) to the perturbation.

where

This is a time correlation function, relating fluctuation of \({\langle O \rangle}\) at time t to the strength of the perturbation \({A}\) at some earlier time \({t'}\). \({K(OA\vert t-t')}\) is defined only for \({t'<t}\), in accordance with the principle of causality. Thus, it is a function only of the difference \({\tau = t - t'}\). Recalling the definitions of the Fourier transform \({f(\omega)}\):

Then instead of 8, we have:

Requiring \({H'(t)}\) to be Hermitian,

Now,

Which, upon combining the expressions for \({H'(t)}\) so as to `Hermitize’ the expression:

Thus

with \({F(t)}\) real. Instead of working in the time domain, we may also consider the response in terms of a single oscillatory perturbation. This means that

To ensure \({H'(\omega)}\) builds smoothly from zero at \({t \rightarrow -\infty}\), we can introduce a convergence factor \({e^{\eta t}}\) with the initial condition \({c_0 = 1}\) and \({c_n = 0}\), which gives:

Then, collecting terms of \({\pm \omega}\):

Finally:

Which is the response function, or frequency-dependent polarizability.