Derivation of the Fluctuation Dissipation Theorem

We previously derived the time evolution operators for both Hamiltonian and Schrodinger equations.

Time evolution operator as a Volterra series

A Heisenberg observable A^H in classical or quantum mechanics will evolve by the time evolution operator

A^H(tb)=𝒰(ta,tb)A^H(ta)

where the derivative is the Liouville operator, and the time evolution operator is the operator exponentiation of the integrated Liouville operator.

dA^H(t)dt=(t)A^H
𝒰(ta,tb)=etatb(t)dt

Now if we substitute back into the derivative

ddt(𝒰(ta,t)A^H(ta))=(t)𝒰(ta,t)A^H(ta)

then integrate both sides

𝒰(ta,tb)A^H(ta)𝒰(ta,ta)A^H(ta)=tatb(t)𝒰(ta,t)A^H(ta)dt

replacing 𝒰(ta,ta) with 1

𝒰(ta,tb)A^H(ta)=A^H(ta)+tatb(t)𝒰(ta,t)A^H(ta)dt

we get the following relation.

𝒰(ta,tb)A^H(ta)=(1+tatb(t)𝒰(ta,t)dt)A^H(ta)

And since it should hold for arbitrary observable, the operators on both sides of the equation must be equal.

𝒰(ta,tb)=1+tatb(t)𝒰(ta,t)dt

Expanding the Volterra series

If we keep substituting into itself, on iteration 0:

𝒰(ta,tb)=1

On iteration 1:

𝒰(ta,tb)=1+tatb(t)dt

Then on iteration 2:

𝒰(ta,tb)=1+tatb(t)dt+tatbtat1(t1)(t2)dt1dt2

Then on iteration 3:

𝒰(ta,tb)=1+tatb(t)dt+tatbtat1(t1)(t2)dt1dt2+tatbtat1tat2(t1)(t2)(t3)dt1dt2dt3

And so on. If we assume it converges the limit should be a solution

𝒰(ta,tb)=1+tatb(t)dt+tatbtat1(t1)(t2)dt1dt2+tatbtat1tat2(t1)(t2)(t3)dt1dt2dt3+

Just for syntax we will group all the higher order terms.

𝒰(ta,tb)=1+tatb(t)dt+𝒪(2)

Time evolution of probability

Classical

If we have a probability distribution ρ(t,q,p) over the phase space we can find the expected value of an observable by

𝔼[A](t)=VA(q,p)ρ(t,q,p)dqdp

And the probability distribution flows according the the continuity equation

dρdt(t,x(t))=x(dxdtρ)

where the state vector x=q1,q2,,p1,p2, so we can split it up and say the below.

dρdt(t,q(t),p(t))=q(dqdtρ)p(dpdtρ)

Use the chain rule

=k(ρqkdqkdt+ρqkdqkdt)k(ρpkdpkdt+ρpkdpkdt)

then substitue in Hamiltons equations.

=k(ρqkHpk+ρ2Hqkpk)k(ρpkHqkρ2Hpkqk)

Cancel those two terms.

=k(ρqkHpkρpkHqk)

It's just the Possion bracket

={H,ρ}

Define an operator that does this

=𝒦(t)ρ(t)

Quantum

If we have a density matrix ρ(t) over the phase space we can find the expected value of an observable by

𝔼[A](t)=Tr(ρ(t)A^)

where the density matrix is defined below.

ρ(t)=kpk|ψk(t)ψk(t)|

If we take the derivative of the density matrix we use the chain rule

dρ(t)dt=kpkddt(|ψk(t)ψk(t)|)=kpk(|dψk(t)dtψk(t)|+|ψk(t)dψk(t)dt|)

and substituting in Schrodingers equation we get

=kpk(iH^(t)|ψk(t)ψk(t)|+|ψk(t)ψk(t)|iH^(t))

and the Hamiltonian is its on Hermetian so

=i(H^(t)ρ(t)=ρ(t)H^(t))

It's just the commutator bracket

=i[H^(t),ρ(t)]

Define an operator that does this

=𝒦(t)ρ(t)

Since the derivative of the probability is

dρdt=𝒦(t)ρ(t)

Make a time evolution operator as the exponentiation

ρ(tb)=etatb𝒦(t)dtρ(ta)=𝒫(ta,tb)ρ(ta)

Interation picture

In the interaction picture we split the Hamiltonian into time independent and a time dependent part.

H=H0+H1(t)

Classical

Define an interaction picture version.

𝒦(t)={H(t),}
={H0(t)+H1(t),}

The Possion bracket distributes over addition.

={H0,}+{H1(t),}

Then define 2 new operators for each term

=𝒦0+𝒦1(t)

Quantum

Define an interaction picture Liouville operator.

𝒦(t)=i[H(t),]
=i[H0+H1(t),]

The commutator bracket distributes over addition.

=i[H0,]i[H1(t),]

Then define 2 new operators for each term

=𝒦0+𝒦1(t)
𝒫0(ta,tb)=e(tbta)𝒦0
𝒫1(ta,tb)=etatb𝒦1(t)dt
𝒦1,I(t)=e(tta)𝒦0𝒦1(t)e(tta)0
𝒫1,I(t)=etatb𝒦1,I(t)dt

and can prove with showing starting condition is the same and derivative is the same

𝒫(ta,tb)=𝒫0(ta,tb)𝒫1,I(ta,tb)

We can do a Volterra series it as

𝒫1,I(ta,tb)=1+tatb𝒦1,I(t)𝒫1,I(ta,t))dt

and expand it

=1+tatb𝒦1,I(t)dt+𝒪(𝒦1,I2)

Then

𝒫(ta,tb)=𝒫0(ta,tb)+tatb𝒫0(ta,tb)𝒦1,I(t)dt+𝒪(𝒦1,I2)
=𝒫0(ta,tb)+tatbe(tbta)𝒦0e(tta)𝒦0𝒦1(t)e(tta)𝒦0dt+𝒪(𝒦1,I2)
=𝒫0(ta,tb)+tatbe(tbt)𝒦0𝒦1(t)e(tta)𝒦0dt+𝒪(𝒦1,I2)
=𝒫0(ta,tb)+tatb𝒫0(t,tb)𝒦1(t)𝒫0(ta,t)dt+𝒪(𝒦1,I2)

Linear response

Consider a closed system where the Hamiltonian is time independent, and is then perturbed by a force along some observable

H(t)=H0h(t)B^

Here

H1(t)=h(t)B^
so

Classical

Here 𝒦1(t) is

𝒦1(t)={h(t)B^,}=h(t){B^,}

So let's define

𝒦B={B^,}

So that

𝒦1(t)=h(t)𝒦B

Quantum

Here 1(t) is

𝒦1(t)=i[h(t)B^,]=h(t)(i[B^,])

So let's define

𝒦B(t)=i[B^(t),]

So that

𝒦1(t)=h(t)𝒦B(t)

So we get

𝒫(ta,tb)=𝒫0(ta,tb)+tatb𝒫0(t,tb)(h(t)𝒦B)𝒫0(ta,t)+𝒪(h2)

And we find that

ρ(tb)ρ0(tb)=tatb𝒫0(t,tb)(h(t)𝒦B)𝒫0(ta,t)ρ(ta)+𝒪(h2)

Expected value of observable at equilibium

Assume that at time ta there was no perturbation and the system is in thermal equilibrium. For all tta the probability distribution is at some equilibrium distribution ρ(t)=ρ0 and there was no force h(t)=0.

Classical

If we have a probability distribution ρ(t,q,p) over the phase space we can find the expected value of an observable by

𝔼[A](t)=VA(q,p)ρ(t,q,p)dqdp

The difference in expected values is

𝔼[A](tb)𝔼0[A](tb)=Vρ(tb)AdqdpVρ0(tb)Adqdp
=V(ρ(tb)ρ0(tb))Adqdp

Then put in the Volterra series

=V(tatb𝒫0(t,tb)(h(t)𝒦B)𝒫0(ta,t)ρ(ta)dt)Adqdp+𝒪(h2)

and move the integrals around

=tatbh(t)V(𝒫0(t,tb)𝒦B𝒫0(ta,t)ρ(ta))Adqdpdt+𝒪(h2)

Call the factor with first order of h the linear response function

χAB(ttb)=V(𝒫0(t,tb)𝒦B𝒫0(ta,t)ρ(ta))Adqdp

put in the equilbrium distribution ρ(ta)=ρ0(q,p)=eH0(q,p)kBTZ

χAB(ttb)=V(𝒫0(t,tb)𝒦B𝒫0(ta,t)ρ0)Adqdp

And by definition the distribution shouldn't evolve under the unberturbed Hamiltonian (you can prove this).

χAB(ttb)=V(𝒫0(t,tb)𝒦Bρ0)Adqdp

so its like

=V(𝒫0(t,tb)iBqiρ0piBpiρ0)qi)Adpdq
=1kBTV𝒫0(t,tb)(eH0(q,p)kBTZ(iBqiH0piBpiH0qi))Adpdq
=1kBTV𝒫0(t,tb)(ρ00B))Adpdq

where 0={,H0} is a Liouville operator for only the equilibrium Hamiltonian. The time evolution operator distributes over multiplication.

=1kBTV(𝒫0(t,tb)ρ0)(𝒫0(t,tb)0B)Adpdq
=1kBTVρ0𝒫0(t,tb)(0B)Adpdq

And since 𝒦0(t)=0(t), and operators commute with their own exponent:

=1kBTVρ0(0e(tbt)0B)Adpdq
=1kBTVρ0(B˙0(ttb))Adpdq

Where B˙ is time shifted to measure ttb in the future. And this is just the expected value at equilibrium

=1kBT𝔼0[B˙(ta+ttb)A(ta)]

So

χAB(ttb)=1kBTddt𝔼0[B(ta+ttb)A(ta)]
χAB(t)=1kBTddt𝔼0[B(ta+t)A(ta)]

Quantum

So the expected value is

𝔼[A](t)=Tr(ρ(t)A^)

And the difference is

𝔼[A](tb)𝔼0[A](tb)=Tr(ρ(tb)A^)Tr(ρ0(tb)A^)
=Tr((ρ(tb)ρ0(tb))A^)

Put the volterra series in there

=Tr(tatb(𝒫0(t,tb)(h(t)𝒦B)𝒫0(ta,t)ρ(ta))A^dt)+𝒪(h2)

Move some stuff around

=tatbh(t)Tr(𝒫0(t,tb)𝒦B(t)𝒫0(ta,t)ρ(ta))A^)dt+𝒪(h2)

Call this first order term the linear response function

χAB(ttb)=Tr((𝒫0(t,tb)𝒦B(t)𝒫0(ta,t)ρ(ta))A^)

Assume equilibrium at ta and put in the equilbium density matrix ρ0=eH^0kBTZ.

=Tr((𝒫0(t,tb)𝒦B(t)𝒫0(ta,t)ρ0)A^)

And you can check that it doesn't change in time from the equilbrium Hamiltonian

=Tr((𝒫0(t,tb)𝒦B(t)ρ0)A^)

Then expand

=Tr(i(𝒫0(t,tb)(ρ0BBρ0))A^)

and because 𝒦(t)=(t)

=Tr(i(ρ0B(ttb)B(ttb)ρ0)A^)

split

=i(Tr(eH^0kBTB(ttb)A^)Tr(B(ttb)eH^0kBTA^))

and remember the cycle property of traces, and cycle the second one

=i(Tr(eH^0kBTB(ttb)A^)Tr(eH^0kBTA^B(ttb)))
=iTr(eH^0kBT[A,B(ttb)])

so

χAB(ttb)=iTr(eH^0kBT[A(ta),B(ta+ttb)])
χAB(t)=iTr(eH^0kBT[A(ta),B(ta+t)])
χAB(t)=i𝔼0[[A(ta),B(ta+t)]]

This is the Kubo formula


The time evolution operator

𝒰(ta,tb)A^=e(tbta)iH^0A^e(tbta)iH^0

Next we consider imaginary time

=i(Tr(eH^0kBTB(t)A^)Tr(eH^0kBTA^B(t)))
=i(Tr(eH^0kBTB(t)eH^0kBTeH^0kBTA^)Tr(eH^0kBTA^B(t)))
=i(Tr(B(tikBT)eH^0kBTA^)Tr(eH^0kBTA^B(t)))

Then cycle the trace again

χAB(t)=i(Tr(eH^0kBTA^B(tikBT))Tr(eH^0kBTA^B(t)))

fourier transform (linear response must be zero for negative times so the integral can start at zero)

χAB(ω)=0eiωti(Tr(eH^0kBTA^B(tikBT))Tr(eH^0kBTA^B(t)))dt
χAB(ω)=i(0eiωtTr(eH^0kBTA^B(tikBT))dt0eiωtTr(eH^0kBTA^B(t))dt)

And we can remove the phase shift in the first term

=i(eωkBT0eiωtTr(eH^0kBTA^B(t))dt0eiωtTr(eH^0kBTA^B(t))dt)
=i(eωkBT1)0eiωtTr(eH^0kBTA^B(t))dt
=i(eωkBT1)0eiωt𝔼0[A^B(t)]dt

Imaginary part

Classical

Now take the fourier transform

χAB(ω)=1kBTeiωtddt𝔼0[B(ta+t)A(ta)]dt

It shouldn't work for negative values so

=1kBT0eiωtddt𝔼0[B(ta+t)A(ta)]dt

integrate by parts

1kBT(eiω𝔼0[B(ta+)A(ta)]e0𝔼0[B(ta)A(ta)])+1kBT0(ddteiωt)𝔼0[B(ta+t)A(ta)]dt

In the limit to infinity they should become uncorrelated

1kBT𝔼0[B(ta)A(ta)])+1kBT0(ddteiωt)𝔼0[B(ta+t)A(ta)]dt
1kBT𝔼0[B(ta)A(ta)])iωkBT0eiωt𝔼0[B(ta+t)A(ta)]dt

taking just the imaginary part

Im[χAB(ω)]=ωkBTRe[0eiωt𝔼0[B(ta+t)A(ta)]dt]

Taking just the real part of anything is

Re[f]=f+f2

So doing that gives

=ω2kBT(0eiωt𝔼0[B(ta+t)A(ta)]dt+0eiωt𝔼0[B(ta+t)A(ta)]dt)
=ω2kBT(0eiωt𝔼0[B(ta+t)A(ta)]dt+0eiωt𝔼0[B(tat)A(ta)]dt)

And if A=B then

=ω2kBT(0eiωt𝔼0[A(ta+t)A(ta)]dt+0eiωt𝔼0[A(tat)A(ta)]dt)

And at equlibrium it's time invariant so you can shift both the offsets in the right integral by t.

=ω2kBT(0eiωt𝔼0[A(ta+t)A(ta)]dt+0eiωt𝔼0[A(ta)A(ta+t)]dt)
=ω2kBTeiωt𝔼0[A(ta+t)A(ta)]dt

So we get a final formulat relating the spectral density to the linear respones function.

eiωt𝔼0[A(ta+t)A(ta)]dt=2kBTωIm[χAA(ω)]

Quantum

Taking just the imaginary part of anything is

Im[f]=ff2i

so

Im[χAB(ω)]=12(eωkBT1)((0eiωt𝔼0[A^B(t)]dt)(0eiωt𝔼0[B(t)A^]dt))

And observables are Hermetian

=12(eωkBT1)((0eiωt𝔼0[A^B(t)]dt)+(0eiωt𝔼0[B(t)A^]dt))

Let A=B

Im[χAA(ω)]=12(eωkBT1)((0eiωt𝔼0[A^A(t)]dt)+(0eiωt𝔼0[A(t)A^]dt))
=12(eωkBT1)((0eiωt𝔼0[A^A(t)]dt)+(0eiωt𝔼0[A(t)A^]dt))

And at equilibrium shifting the time shouldn't matter 𝔼0[A(tat)A(ta)]=𝔼0[A(ta)A(ta+t)]

=12(eωkBT1)((0eiωt𝔼0[A^A(t)]dt)+(0eiωt𝔼0[AA^(t)]dt))

Then combined makes

12(eωkBT1)eiωt𝔼0[A^A(t)]dt

And we get

eiωt𝔼0[A^A(t)]dt=2eωkBT1Im[χA(ω)]

Classical limit

If we negative the frequency

S(ω)=S(ω)=eiωt𝔼0[A^(t)A^(t0)]dt=2(1eβω1+1)Im[χ(ω)]

And since χ is a real function the Imaginary component is an odd function

=2(1eβω1+1)Im[χ(ω)]
=2(eβωeβω11)Im[χ(ω)]
=21eβω1Im[χ(ω)]

So in the quantum version unlike the classical version S(ω)S(ω), the spectral density isn't the same when frequency is negated. To fix this we can average the two.

S(ω)+S(ω)2=2(1eβω1+12)Im[χ(ω)]
=(2eβω1+1)Im[χ(ω)]
=eβω+1eβω1Im[χ(ω)]
=eβω2+eβω2eβω2eβω2Im[χ(ω)]

And thats just the formula for hyperbolic cotangent

=coth(βω2)Im[χ(ω)]

for large temepratrues

lim0coth(βω2)=2βω

so

lim0S(ω)+S(ω)2=2βωIm[χ(ω)]

And this recovers the classical version in the classical limit

lim0S(ω)+S(ω)2=2kBTωIm[χ(ω)]