Problem 4-1: Fermionic double-well

Hi guys. Here is a notebook with a bunch of stuff about problem 4.1 
I don't really care what you do with this problem, I was just hoping you explore this notebook a bit to try to understand the differences between the canonical and grand-canonical ensemble. 

I computed some thermodynamic quantities here, but there are many others to compute. So feel free to explore and compute interesting things.

About Fermionic operators in Mathematica

Qulib has a function LoadSpinChainMatrices which loads spin 1/2 operators for a quantum spin chain

LoadSpinChainMatrices[2,"a"];

LoadedMatrices

Using the Jordan-Wigner transformation, this function also loads fermionic operators c[i] and cd[i] = c[i]†

The matrix elements of these operators can look a bit weird.
But all that matters is that they satisfy the fermionic algebra.

{(Subscript[c, iσ]),(Subscript[c, jσ']^†)} = (Subscript[δ, ij])(Subscript[δ, σσ'])

c[1]//mf;
c[2]//mf;

c[1]\[ConjugateTranspose].c[1]+c[1].c[1]\[ConjugateTranspose]//mf;

c[1]\[ConjugateTranspose].c[2]+c[2].c[1]\[ConjugateTranspose]//mf;

As a sanity check, we can also use this to verify that the Bogoliubov transformation we played with in class does lead to fermionic operators.

These are the formulas we wrote in class

b[1] = Cos[(ϕ/2)] c[1]-I Sin[(ϕ/2)] c[2]\[ConjugateTranspose];
b[2] = I Sin[(ϕ/2)] c[1]\[ConjugateTranspose]+Cos[(ϕ/2)] c[2];

The matrix elements are again unintuitive:

b[1]//mf;

But all that matters is the algebra:

b[1].b[1]\[ConjugateTranspose]+b[1]\[ConjugateTranspose].b[1]//cf//mf;
b[2].b[2]\[ConjugateTranspose]+b[2]\[ConjugateTranspose].b[2]//cf//mf;

b[1].b[2]\[ConjugateTranspose]+b[2]\[ConjugateTranspose].b[1]//cf//mf;

Writing down the Hamiltonian

This notebook uses qulib to generate 4 fermionic matrices labeled by (i,σ), where i ∈ {L,R} and σ ∈ {1,-1}

LoadSpinChainMatrices[4,"a"];
ivals = {"L","R"};
σvals = {"+","-"};

c["L","+"] = c[1];
c["L","-"] = c[2];
c["R","+"] = c[3];
c["R","-"] = c[4];

(* This defines the number operator (Subscript[n, iσ]) *)
Do[n[i,σ] = c[i,σ]\[ConjugateTranspose].c[i,σ],{i,ivals},{σ,σvals}]

(* This is the total number of fermions *)
𝒩 = Sum[n[i,σ],{i,ivals},{σ,σvals}];

Now this is how one may write the Hamiltonian in Problem 4.1.
Here and henceforth I will choose J = 1

Clear[J,U]
J = 1;
H = U Sum[n[i,"+"].n[i,"-"],{i,ivals}] - J Sum[c["L",σ]\[ConjugateTranspose].c["R",σ] +c["R",σ]\[ConjugateTranspose].c["L",σ],{σ,σvals}]//mf;

The matrix elements are again confusing. 
But we can compute stuff such as eigenvalues and so on.

eigs=Eigenvalues[H]

vecs=FullSimplify[Normalize/@Eigenvectors[H],U>0];

Grand-canonical ensemble

As we discussed in class, the thermodynamics of the grand-canonical ensemble will coincide with that of the canonical ensemble in the thermodynamic limit. 

Here this is not the case: since we only have 2 sites and 2 spin values there can be at most 4 particles present.

Thus, the grand-canonical and canonical ensembles in this case give different physics.

The grand-canonical ensemble  should thus be interpreted as a particle reservoir  which allows for (OverHat[N]) to fluctuate, such that only ⟨(OverHat[N])⟩ = N is fixed.

About the chemical potential

In order to allow for the number of particles may fluctuate,  we need is to include a chemical potential

ρ = MatrixExp[-(H - μ 𝒩)/T];
Z = Tr[ρ];
ρ = ρ/Z;

The average number of particles is now a funciton of μ

Nave = Tr[𝒩.ρ]//Simplify

Here is ⟨(OverHat[N])⟩ vs. μ

Plot[Nave/.{U->3,T->0.1},{μ,-7,7},GridLines->{None,{0,1,2,3,4}},
FrameLabel->{"μ","⟨N⟩"}
]

But what we want is the other way around: we want μ as a function of ⟨(OverHat[N])⟩ = N

So here is a function which finds the inverse μ(U,T,N)

(* μ0 is an initial guess for the root-finding function *)
Clear[FindChemicalPotential];
FindChemicalPotential[U_,T_,ℕ_,μ0_:0]:=Module[{nave = (2 (2+(E^((3 U+(Sqrt[16+(U^2)])-4 μ)/(2 T)))+(E^((-1+2 U-3 μ)/T))+(E^((U-2 μ)/T))+3 (E^((-1+U-μ)/T))+3 (E^((2 (U-μ))/T))+3 (E^((1+U-μ)/T))+(E^(-((-3 U+(Sqrt[16+(U^2)])+4 μ)/(2 T))))+(E^((1/T)+((2 U)/T)-((3 μ)/T)))))/(1+(E^((3 U+(Sqrt[16+(U^2)])-4 μ)/(2 T)))+2 (E^((-1+2 U-3 μ)/T))+(E^((U-2 μ)/T))+(E^((2 (U-2 μ))/T))+2 (E^((-1+U-μ)/T))+3 (E^((2 (U-μ))/T))+2 (E^((1+U-μ)/T))+(E^(-((-3 U+(Sqrt[16+(U^2)])+4 μ)/(2 T))))+2 (E^((1/T)+((2 U)/T)-((3 μ)/T))))},μ/.FindRoot[nave==ℕ,{μ,μ0}]]

FindChemicalPotential[1.9,0.5,2]

Here is the half-filling case ⟨(OverHat[N])⟩ = 2. 
I plot μ vs. U for different T

tab = Table[{U,FindChemicalPotential[U,T,2]},{T,{0.05,0.1,0.5}},{U,0,3,0.3}];
ListLinePlot[tab,
PlotLegends->LineLegend[Automatic,{"T = 0.05","T = 0.1","T = 0.5"}],
FrameLabel->{"U","μ"},
PlotLabel->"⟨(OverHat[N])⟩ = 2"
]

This numerics suggests that when ⟨(OverHat[N])⟩ = 2, we should have μ = U/2 for any temperature.
This is the half-filling case (one particle per site)

Nave/.μ->U/2//Simplify

Voilá

For other values of N this is no longer true and μ becomes a non-linear function of U.
Here is ⟨(OverHat[N])⟩ = 1

tab = Table[{U,FindChemicalPotential[U,T,1]},{T,{0.3}},{U,0,5,0.3}];
ListLinePlot[tab,
PlotLegends->LineLegend[Automatic,{"T = 0.3"}],
FrameLabel->{"U","μ"},
PlotLabel->"⟨(OverHat[N])⟩ = 1"
]

Technical point:

When N = 3, the function FindChemicalPotential becomes buggy at low temperatures.
So there is an option 4th element  in "FindChemicalPotential" which allows you to specify an initial guess for the root finding.

FindChemicalPotential[1,0.1,3]

FindChemicalPotential[1,0.1,3,3]

tab = Quiet@Table[{U,FindChemicalPotential[U,T,3]},{T,{0.1,0.3}},{U,0,10,0.3}];
ListLinePlot[tab]

tab = Quiet@Table[{U,FindChemicalPotential[U,T,3,2U]},{T,{0.05,0.1,0.3}},{U,0,10,0.3}];
ListLinePlot[tab]

Thermodynamic quantities

Now that we have the chemical potential, we can compute other quantities as usual.
All we need to remember is that the chemical potential should always be replaced by the function μ(U,T,N)

To illustrate this, I will compute here only the internal energy:

Have = Tr[H.ρ]//Simplify

In this formula μ is still left implicit. To make plots of ⟨H⟩ we substitute μ with FindChemicalPotential

ℕ = 1;
tab = Table[{T,Have/.μ->FindChemicalPotential[U,T,ℕ]},{U,{0.1,1.0,3.0}},{T,0.05,0.7,0.05}];
ListLinePlot[tab,FrameLabel->{"T","⟨H⟩"},
PlotLabel->"⟨(OverHat[N])⟩ = 1",PlotRange->All,
PlotLegends->LineLegend[Automatic,{"U = 0.1", "U = 1", "U = 3"}]]

ℕ = 2;
tab = Table[{T,Have/.μ->FindChemicalPotential[U,T,ℕ]},{U,{0.1,1.0,3.0}},{T,0.01,1,0.05}];
ListLinePlot[tab,FrameLabel->{"T","⟨H⟩"},
PlotLabel->"⟨(OverHat[N])⟩ = 2",
PlotLegends->LineLegend[Automatic,{"U = 0.1", "U = 1", "U = 3"}]]

ℕ = 3;
tab = Table[{T,Have/.μ->FindChemicalPotential[U,T,ℕ,2U]},{U,{0.1,1.0,3.0}},{T,0.05,1,0.01}];
ListLinePlot[tab,FrameLabel->{"T","⟨H⟩"},
PlotLabel->"⟨(OverHat[N])⟩ = 3",PlotRange->All,
PlotLegends->LineLegend[Automatic,{"U = 0.1", "U = 1", "U = 3"}]]

ℕ = 4;
tab = Table[{T,Have/.μ->FindChemicalPotential[U,T,ℕ,2U]},{U,{0.1,1.0,3.0}},{T,0.05,1,0.01}];
ListLinePlot[tab,FrameLabel->{"T","⟨H⟩"},
PlotLabel->"⟨(OverHat[N])⟩ = 4",PlotRange->All,
PlotLegends->LineLegend[Automatic,{"U = 0.1", "U = 1", "U = 3"}]]

You can now proceed similarly and compute other thermodynamic quantities

Canonical ensemble

The Hamiltonian H conserves the number of fermions, so we can order H  into blocks with  well defined number of particles.

Here is now the number operator (OverHat[N]) is structured (this depends on the way we ordered our fermionic operators)

𝒩//mf;

We can define an ordering operation that groups states with the same number of particles. We can then define a permutation of the columns and rows of H so as to order it in blocks with well defined N.

ord=Ordering@Diagonal@𝒩
Hord = H[[ord,ord]];

The original H is spread out:

MatrixPlot[H]

But (Subscript[H, ord]) is now block diagonal into sub-spaces of well-defined numbers of particles.

MatrixPlot[Hord,Mesh->{{1,5,11,15},{1,5,11,15}},PlotRangePadding->None]

We can also define sub-Hamiltonians living in each sub-space.
To do that, we define projectors into sub-spaces of well defined particle number.
These are called isometries.

ℕs = Tally[Diagonal@𝒩];
Clear[Π];
Do[
pos = Flatten@Position[Diagonal@𝒩,ℕ[[1]]];
Π[ℕ[[1]]] = Sum[out[Eye[ℕ[[2]]][[i]],Eye[16][[pos[[i]]]]],{i,Length@pos}];
,{ℕ,ℕs}]
(*Clear[pos];*)

For instance Π[3] projects into the sub-space with N = 3 particles

Π[0]//mf;
Π[1]//mf;
Π[2]//mf;
Π[3]//mf;
Π[4]//mf;

We can now defined reduced Hamiltonians living on each sub-space

Clear[Hred]
Do[Hred[ℕ] = Π[ℕ].H.Π[ℕ]\[ConjugateTranspose],{ℕ,0,4}]

Hred[0]//mf;
Hred[1]//mf;
Hred[2]//mf;
Hred[3]//mf;
Hred[4]//mf;

Note how at N = 1, (Subscript[H, 1]) does not depend on U, because there are no interactions.

If the number of particles is fixed, we can then study the thermodynamics properties at each sub-space.

Eigenvalues[Hred[0]]

Eigenvalues[Hred[1]]

Eigenvalues[Hred[2]]

Eigenvalues[Hred[3]]

Eigenvalues[Hred[4]]

Here is the canonical density matrix ρ = (e^(-βH)) at each sub-space

Do[
ρcan[ℕ] = MatrixExp[- Hred[ℕ]/T];
Zcan[ℕ] = Tr[ρcan[ℕ]]//Simplify;
ρcan[ℕ] = ρcan[ℕ]/Zcan[ℕ];
,{ℕ,0,4}]

Everything we know from thermodynamics will now hold on each sub-space independently:
For instance, here is the internal energy

Do[
HaveCan[ℕ] = Tr[Hred[ℕ].ρcan[ℕ]]//Simplify,
{ℕ,0,4}]

HaveCan[0]
HaveCan[1]
HaveCan[2]
HaveCan[3]
HaveCan[4]

Note how the average energy is 2U in the N = 4 subspace. 
In this subspace all available sites are filled, so thermal fluctuations don't do anything.

Table[Plot[HaveCan[ℕ]/.U->1,{T,0.01,0.7},ImageSize->200],{ℕ,0,4}]

We can compare this with the grand-canonical ensemble

ℕ = 0;
tab = Table[{T,Have/.μ->FindChemicalPotential[U,T,ℕ]},{U,{0.3}},{T,0.01,1,0.05}];
Show[
Plot[HaveCan[ℕ]/.{β->1/T,U->0.3},{T,0.01,1}],
ListLinePlot[tab,PlotMarkers->Automatic],
PlotLabel->"⟨(OverHat[N])⟩ = "<>ToString@ℕ,
FrameLabel->{"T","⟨H⟩"}
]

ℕ = 1;
tab = Table[{T,Have/.μ->FindChemicalPotential[U,T,ℕ]},{U,{0.3}},{T,0.01,1,0.05}];
Show[
Plot[HaveCan[ℕ]/.{β->1/T,U->0.3},{T,0.01,1}],
ListLinePlot[tab,PlotMarkers->Automatic],
PlotLabel->"⟨(OverHat[N])⟩ = "<>ToString@ℕ,
FrameLabel->{"T","⟨H⟩"}
]

ℕ = 2;
tab = Table[{T,Have/.μ->FindChemicalPotential[U,T,ℕ]},{U,{0.3}},{T,0.01,1,0.05}];
Show[
Plot[HaveCan[ℕ]/.{β->1/T,U->0.3},{T,0.01,1}],
ListLinePlot[tab,PlotMarkers->Automatic],
PlotLabel->"⟨(OverHat[N])⟩ = "<>ToString@ℕ,
FrameLabel->{"T","⟨H⟩"}
]

ℕ = 3;
tab = Table[{T,Have/.μ->FindChemicalPotential[U,T,ℕ,2 0.3]},{U,{0.3}},{T,0.01,1,0.05}];
Show[
Plot[HaveCan[ℕ]/.{β->1/T,U->0.3},{T,0.01,1}],
ListLinePlot[tab,PlotMarkers->Automatic],
PlotLabel->"⟨(OverHat[N])⟩ = "<>ToString@ℕ,
FrameLabel->{"T","⟨H⟩"}
]

ℕ = 4;
tab = Table[{T,Have/.μ->FindChemicalPotential[U,T,ℕ,2 0.3]},{U,{0.3}},{T,0.01,1,0.05}];
Show[
Plot[HaveCan[ℕ]/.{β->1/T,U->0.3},{T,0.01,1}],
ListLinePlot[tab,PlotMarkers->Automatic],
PlotLabel->"⟨(OverHat[N])⟩ = "<>ToString@ℕ,
FrameLabel->{"T","⟨H⟩"}
]

They don't coincide. But they don't have to. 
The two ensembles only coincide in the thermodynamic limit.