Fluctuation Dissipation Theorem

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

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)=e∫tatbℒ(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+∫tatb∫tat1ℒ(t1)ℒ(t2)dt1dt2

Then on iteration 3:

𝒰(ta,tb)=1+∫tatbℒ(t)dt+∫tatb∫tat1ℒ(t1)ℒ(t2)dt1dt2+∫tatb∫tat1∫tat2ℒ(t1)ℒ(t2)ℒ(t3)dt1dt2dt3

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

𝒰(ta,tb)=1+∫tatbℒ(t)dt+∫tatb∫tat1ℒ(t1)ℒ(t2)dt1dt2+∫tatb∫tat1∫tat2ℒ(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(Γ→)ρ(t,Γ→)dΓ→

where the state vector Γ→=⟨q→,p→⟩=⟨q1,q2,…,p1,p2,…⟩. The probability distribution flows in phase space according the the continuity equation.

dρdt(t,Γ→(t))=−∇Γ⋅(dΓ→dtρ)
dρdt(t,Γ→(t))=−∇q⋅(dq→dtρ)−∇p⋅(dp→dtρ)

Use the chain rule.

=−∑k(∂ρ∂qkdqkdt+ρ∂∂qkdqkdt)−∑k(∂ρ∂pkdpkdt+ρ∂∂pkdpkdt)

Then substitue the time evolution of the states with Hamiltons equations.

=−∑k(∂ρ∂qk∂H∂pk+ρ∂2H∂qk∂pk)−∑k(−∂ρ∂pk∂H∂qk−ρ∂2H∂pk∂qk)

The terms involving second derivatives of the Hamiltonian can be canceled.

=−∑k(∂ρ∂qk∂H∂pk−∂ρ∂pk∂H∂qk)

The right hand side is just the Possion bracket.

={H,ρ}

We can define an operator on the probability distribution that does this.

𝒦(t)={H(t),⋅}

Then the probility distribution evolves by the operator.

dρdt=𝒦(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)|

Taking 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|)

We can then substitute Schrodingers equation for the derivatives.

=∑kpk(−iℏH^(t)|ψk(t)⟩⟨ψk(t)|+|ψk(t)⟩⟨ψk(t)|iℏH^†(t))

The Hamiltonian is its own Hermetian. Then simplify.

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

The right hand side is just the commutator bracket.

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

We can define an operator on the density matrix that does this.

𝒦(t)=−iℏ[H^(t),⋅]

Then the density matrix evolves by the operator..

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

Since the derivative of the probability is

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

a time evolution operator can be defined as the exponentiation.

ρ(tb)=e∫tatb𝒦(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)

Then we can define a probability time evolution for each part of the Hamiltonian.

𝒫0(ta,tb)=e(tb−ta)𝒦0
𝒫1(ta,tb)=e∫tatb𝒦1(t′)dt′

Then define an interaction picture version of the time dependent part.

𝒦1,I(t)=e−(t−ta)𝒦0𝒦1(t)e(t−ta)𝒦0
𝒫1,I(t)=e∫tatb𝒦1,I(t′)dt′

One can prove by showing the starting condition is the same and derivatives are the same that the below is true.

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

We can set up a Volterra series for the operator as

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

and expand it as before.

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

Then multiplying by the time independent part we get the exact probability time evolution operator.

𝒫(ta,tb)=𝒫0(ta,tb)+∫tatb𝒫0(ta,tb)𝒦1,I(t)dt+𝒪(𝒦1,I2)
=𝒫0(ta,tb)+∫tatbe(tb−ta)𝒦0e−(t−ta)𝒦0𝒦1(t)e(t−ta)𝒦0dt+𝒪(𝒦1,I2)
=𝒫0(ta,tb)+∫tatbe(tb−t)𝒦0𝒦1(t)e(t−ta)𝒦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)=H0−h(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 t≤ta 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→)dq→dp→

The difference in expected values is

𝔼[A](tb)−𝔼0[A](tb)=∫Vρ(tb)Adq→dp→−∫Vρ0(tb)Adq→dp→
=∫V(ρ(tb)−ρ0(tb))Adq→dp→

Then put in the Volterra series

=∫V(∫tatb𝒫0(t,tb)(−h(t)𝒦B)𝒫0(ta,t)ρ(ta)dt)Adq→dp→+𝒪(h2)

and move the integrals around

=∫tatbh(t)∫V−(𝒫0(t,tb)𝒦B𝒫0(ta,t)ρ(ta))Adq→dp→dt+𝒪(h2)

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

χAB(t−tb)=−∫V(𝒫0(t,tb)𝒦B𝒫0(ta,t)ρ(ta))Adq→dp→

put in the equilbrium distribution ρ(ta)=ρ0(q→,p→)=e−H0(q→,p→)kBTZ

χAB(t−tb)=−∫V(𝒫0(t,tb)𝒦B𝒫0(ta,t)ρ0)Adq→dp→

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

χAB(t−tb)=−∫V(𝒫0(t,tb)𝒦Bρ0)Adq→dp→

so its like

=−∫V(𝒫0(t,tb)∑i∂B∂qi∂ρ0∂pi−∂B∂pi∂ρ0)∂qi)Adp→dq→
=1kBT∫V𝒫0(t,tb)(e−H0(q→,p→)kBTZ(∑i∂B∂qi∂H0∂pi−∂B∂pi∂H0∂qi))Adp→dq→
=1kBT∫V𝒫0(t,tb)(ρ0ℒ0B))Adp→dq→

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

=1kBT∫V(𝒫0(t,tb)ρ0)(𝒫0(t,tb)ℒ0B)Adp→dq→
=1kBT∫Vρ0𝒫0(t,tb)(ℒ0B)Adp→dq→

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

=1kBT∫Vρ0(ℒ0e−(tb−t)ℒ0B)Adp→dq→
=1kBT∫Vρ0(B˙0(t−tb))Adp→dq→

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

=1kBT𝔼0[B˙(ta+t−tb)A(ta)]

So

χAB(t−tb)=1kBTddt𝔼0[B(ta+t−tb)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(t−tb)=−Tr((𝒫0(t,tb)𝒦B(t)𝒫0(ta,t)ρ(ta))A^)

Assume equilibrium at ta and put in the equilbium density matrix ρ0=e−H^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)(ρ0B−Bρ0))A^)

and because 𝒦(t)=−ℒ(t)

=Tr(iℏ(ρ0B(t−tb)−B(t−tb)ρ0)A^)

split

=iℏ(Tr(e−H^0kBTB(t−tb)A^)−Tr(B(t−tb)e−H^0kBTA^))

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

=iℏ(Tr(e−H^0kBTB(t−tb)A^)−Tr(e−H^0kBTA^B(t−tb)))
=−iℏTr(e−H^0kBT[A,B(t−tb)])

so

χAB(t−tb)=−iℏTr(e−H^0kBT[A(ta),B(ta+t−tb)])
χAB(t)=−iℏTr(e−H^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(tb−ta)iℏH^0A^e(tb−ta)−iℏH^0

Next we consider imaginary time

=iℏ(Tr(e−H^0kBTB(t)A^)−Tr(e−H^0kBTA^B(t)))
=iℏ(Tr(e−H^0kBTB(t)eH^0kBTe−H^0kBTA^)−Tr(e−H^0kBTA^B(t)))
=iℏ(Tr(B(t−iℏkBT)e−H^0kBTA^)−Tr(e−H^0kBTA^B(t)))

Then cycle the trace again

χAB(t)=iℏ(Tr(e−H^0kBTA^B(t−iℏkBT))−Tr(e−H^0kBTA^B(t)))

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

χAB(ω)=∫0∞e−iωtiℏ(Tr(e−H^0kBTA^B(t−iℏkBT))−Tr(e−H^0kBTA^B(t)))dt
χAB(ω)=iℏ(∫0∞e−iωtTr(e−H^0kBTA^B(t−iℏkBT))dt−∫0∞e−iωtTr(e−H^0kBTA^B(t))dt)

And we can remove the phase shift in the first term

=iℏ(eℏωkBT∫0∞e−iωtTr(e−H^0kBTA^B(t))dt−∫0∞e−iωtTr(e−H^0kBTA^B(t))dt)
=iℏ(eℏωkBT−1)∫0∞e−iωtTr(e−H^0kBTA^B(t))dt
=iℏ(eℏωkBT−1)∫0∞e−iωt𝔼0[A^B(t)]dt

Imaginary part

Classical

Now take the fourier transform

χAB(ω)=1kBT∫−∞∞e−iωtddt𝔼0[B(ta+t)A(ta)]dt

It shouldn't work for negative values so

=1kBT∫0∞e−iωtddt𝔼0[B(ta+t)A(ta)]dt

integrate by parts

1kBT(e−iω∞𝔼0[B(ta+∞)A(ta)]−e0𝔼0[B(ta)A(ta)])+1kBT∫0∞(ddte−iωt)𝔼0[B(ta+t)A(ta)]dt

In the limit to infinity they should become uncorrelated

−1kBT𝔼0[B(ta)A(ta)])+1kBT∫0∞(ddte−iωt)𝔼0[B(ta+t)A(ta)]dt
−1kBT𝔼0[B(ta)A(ta)])−iωkBT∫0∞e−iωt𝔼0[B(ta+t)A(ta)]dt

taking just the imaginary part

Im[χAB(ω)]=−ωkBTRe[∫0∞e−iωt𝔼0[B(ta+t)A(ta)]dt]

Taking just the real part of anything is

Re[f]=f+f∗2

So doing that gives

=−ω2kBT(∫0∞e−iωt𝔼0[B(ta+t)A(ta)]dt+∫0∞eiωt𝔼0[B(ta+t)A(ta)]dt)
=−ω2kBT(∫0∞e−iωt𝔼0[B(ta+t)A(ta)]dt+∫−∞0e−iωt𝔼0[B(ta−t)A(ta)]dt)

And if A=B then

=−ω2kBT(∫0∞e−iωt𝔼0[A(ta+t)A(ta)]dt+∫−∞0e−iωt𝔼0[A(ta−t)A(ta)]dt)

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

=−ω2kBT(∫0∞eiωt𝔼0[A(ta+t)A(ta)]dt+∫−∞0eiωt𝔼0[A(ta)A(ta+t)]dt)
=−ω2kBT∫−∞∞eiω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]=f−f∗2i

so

Im[χAB(ω)]=12ℏ(eℏωkBT−1)((∫0∞e−iωt𝔼0[A^B(t)]dt)−(−∫0∞eiωt𝔼0[B†(t)A^†]dt))

And observables are Hermetian

=12ℏ(eℏωkBT−1)((∫0∞e−iωt𝔼0[A^B(t)]dt)+(∫0∞eiωt𝔼0[B(t)A^]dt))

Let A=B

Im[χAA(ω)]=12ℏ(eℏωkBT−1)((∫0∞e−iωt𝔼0[A^A(t)]dt)+(∫0∞eiωt𝔼0[A(t)A^]dt))
=12ℏ(eℏωkBT−1)((∫0∞e−iωt𝔼0[A^A(t)]dt)+(∫−∞0e−iωt𝔼0[A(−t)A^]dt))

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

=12ℏ(eℏωkBT−1)((∫0∞e−iωt𝔼0[A^A(t)]dt)+(∫−∞0e−iωt𝔼0[AA^(t)]dt))

Then combined makes

12ℏ(eℏωkBT−1)∫−∞∞e−iωt𝔼0[A^A(t)]dt

And we get

∫−∞∞e−iωt𝔼0[A^A(t)]dt=2ℏeℏωkBT−1Im[χA(ω)]

Classical limit

If we negative the frequency

S′(ω)=S(−ω)=∫−∞∞e−iω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βℏω−1−1)Im[χ(ω)]
=2ℏ1eβℏω−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βℏω2−e−βℏω2Im[χ(ω)]

And thats just the formula for hyperbolic cotangent

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

for large temepratrues

limℏ→0⁡ℏcoth⁡(βℏω2)=2βω

so

limℏ→0⁡S(ω)+S(−ω)2=2βωIm[χ(ω)]

And this recovers the classical version in the classical limit

limℏ→0⁡S(ω)+S(−ω)2=2kBTωIm[χ(ω)]