Load MathQ - beta

SetDirectory[NotebookDirectory[]];

nb=CreateDocument[Get["MathQ-beta.nb"], WindowSelected->False];
NotebookEvaluate@nb;
NotebookClose[nb];

Lesson 1 - Building a Lattice

Basic lattices

Built-in lattice presets

We have a number of basic lattice presets (in 1D, 2D and 3D), and some more elaborate lattice models (we will not discuss these)

LatticePresets

lat = BuildLattice["Honeycomb"]

List@@lat

lat["Sites"]

lat["Tags"]

lat["LatticeDimensions"]

lat["BravaisVectors"]

A useful command is PlotLattice. By default it plots the unitcell and (in dimmer color) the neighboring cells along the Bravais vectors (a1, a2, and a1+a2 here).
Note that PlotLattice accepts any option from Graphics (see Options[PlotLattice] for complete list)

PlotLattice[lat]
PlotLattice[lat, Axes->True,Frame->True,BaseStyle->FontSize->15]

We can also choose to plot more (or less) unit cells

PlotLattice[lat, UnitCells->10, Axes->True, Frame->True,BaseStyle->FontSize->15, BoxCenterAtOrigin->True]

The lattice constant is always 1 in the built-in presets. We can change that (in any lattice) like this

lat=BuildLattice["Honeycomb",LatticeConstant->2]
lat["BravaisVectors"]
Norm/@%
PlotLattice[lat, Axes->True,Axes->True, Frame->True,BaseStyle->FontSize->15]

Custom lattices

We can provide arbitrary bravais vectors and basis sites for any lattice

With[{
bravais={{1,0,0},{0,1,0},{0,0,1}}, sublats = {"A"->{{0,0,0}},"B"->{{0,.5,0},{.6,0,.6}}}},
lat = BuildLattice[{bravais,sublats}]
]

PlotLattice[lat, Axes->True, Boxed->True]

PlotLattice[lat, Axes->True, Boxed->True, Labels->"Numbers"]

PlotLattice[lat, Axes->True, Boxed->True, Labels->"Tags"]

Derived lattices

Superlattices

lat = BuildLattice["Square"]

A supercell that is 10 times larger in all directions

lat2=BuildLattice[lat, Supercell->10]

PlotLattice[lat2]

A supercell that is 3 times larger along the first bravais vector, and 5 times larger along the second

lat3=BuildLattice[lat, Supercell->{3,5}]

PlotLattice[lat3]

A supercell that is generated by a1’ = 4a1 - 4a2 and a2’ = 5a2

lat4=BuildLattice[lat, Supercell->{{4,-4},{0,5}}]

PlotLattice[lat4]

We can also do it in one go

PlotLattice[BuildLattice["Square", Supercell->{{4,-4},{0,5}}]]

Regions

lat = BuildLattice["Square"]

PlotLattice[lat]

Note: We may draw cells boxes with their center (instead of their vertex) at the origin. We can also specify a line style

PlotLattice[lat, BoxCenterAtOrigin->True, CellBoxStyle->Darker[Green]]

We now wish to “fill” a bounded region with this periodic lattice. The result is a non-periodic (bounded) lattice

myregion=Function[{r},Norm[r]<30];
lat2=BuildLattice[lat, Bounded->{Region->myregion}]
PlotLattice[lat2]

A more complicated example

lat = BuildLattice["BCC", Supercell->{3,1,1}]

PlotLattice[lat, UnitCells->5]

Note: In 3D we can render a “pretty” lattice using a finite SiteRadius (caution: unstable for large lattices)

PlotLattice[lat, UnitCells->5, SiteRadius->.4, BoxedCells->None]

We will make lat bounded along two axes (second and third), leaving the first one unbounded

myregion = Function[{r},Norm[r[[2;;3]]]<8];
lat2=BuildLattice[lat, Bounded->{Axes->{2,3}, Region->myregion}]
PlotLattice[lat2, Boxed->True,SiteRadius->.4]

Reshaping unit cells

myregion = Function[{r},Norm[r[[2;;3]]]<8];
lat2=BuildLattice["BCC", Supercell->{4,1,1},Bounded->{Axes->{2,3}, Region->myregion}]
PlotLattice[lat2, Boxed->True,SiteRadius->.4]

We can choose the unit cell of a periodic lattice in many ways that generate, upon translation along the bravais vectors, the same sublattice. We can thus “reshape” the unit cell freely to our preferred shape.

myregion=Function[{r},r[[1]]>0];
lat3=BuildLattice[lat2,Reshape->{Axis->1,Region->myregion}]
PlotLattice[lat3,Boxed->True,SiteRadius->.4]

Deforming lattices

myregion=Function[{r},Norm[r]<10];
lat=BuildLattice["Honeycomb",EmbeddingDimensions->3,Bounded->{Region->myregion}]
PlotLattice[lat, Boxed->True,SiteRadius->.3, SiteColors->{White,Darker[Green]}, Axes->True]

mydeformation=Function[{r},r+4Exp[-(((Norm[r])^2)/(2*(3^2)))]{0,0,1}];
lat2=BuildLattice[lat, Transform->mydeformation]
PlotLattice[lat2, SiteRadius->.3,SiteColors->{White, Darker[Green]}, Boxed->True]

We can pass the site tag as a second argument to spatial functions

mydeformation=Function[{r, tag},r+4Exp[-(((Norm[r])^2)/(2*(3^2)))]{0,0,1}If[tag=="A",1,-1]];
lat2=BuildLattice[lat, Transform->mydeformation]
PlotLattice[lat2, SiteRadius->.3,SiteColors->{White, Darker[Green]}, Boxed->True]

Deformations of periodic lattices are applied only to the base unit cell, and are therefore extended periodically

lat=BuildLattice["Honeycomb",EmbeddingDimensions->3,Supercell->20]
mydeformation=Function[{r},r+4Exp[-(((Norm[r-{0,20((Sqrt[3])/2),0}])^2)/(2*(3^2)))]{0,0,1}];
lat2=BuildLattice[lat, Transform->mydeformation]
PlotLattice[lat2, SiteRadius->.3,SiteColors->{White, Darker[Green]}, BoxedCells->False]

Combining lattices

botlat=BuildLattice["Honeycomb",EmbeddingDimensions->3,Tags->{"Ab", "Bb"},
Bounded->{Region->Function[{r},Norm[r]<15]}, 
Transform->Function[{r}, r-{0,0,.5}]]
toplat=BuildLattice["Honeycomb",EmbeddingDimensions->3,Tags->{"At", "Bt"},
Bounded->{Region->Function[{r},Norm[r]<15]}, 
Transform->Function[{r}, (({{Cos[.2], -Sin[.2], 0}, {Sin[.2], Cos[.2], 0}, {0, 0, 1}})).r+{0,0,.5}]]
PlotLattice[{botlat,toplat},SiteRadius->.2, SiteColors->{{Black,Gray},{Purple,White}}, ViewPoint->Top]
PlotLattice[{botlat,toplat},SiteRadius->.2, SiteColors->{{Black,Gray},{Purple,White}}, ViewPoint->Front]

Lattices can be combined into one if they have the same periodicity (same Bravais vectors - none in the example above)

lat=CombineLattices[botlat, toplat]
PlotLattice[lat,SiteRadius->.2, ViewPoint->Top]

An example with shared periodicity

lat1=BuildLattice["Square",Supercell->4]
lat2=BuildLattice["Square",Supercell->2, Transform->Function[{r},2r+{1/2,1/2}]]
lat3=CombineLattices[lat1,lat2]

PlotLattice[{lat1,lat2}]
lat1["BravaisVectors"]
lat2["BravaisVectors"]

lat3=CombineLattices[lat1,lat2]
PlotLattice[lat3]

Attempting to combine lattices with different Bravais vectors throws an “Incompatible lattices” error

lat1=BuildLattice["Square",Supercell->4]
lat2=BuildLattice["Square",Supercell->2]
lat1["BravaisVectors"]
lat2["BravaisVectors"]
lat3=CombineLattices[lat1,lat2]

Lesson 2 - Building a System

One orbital per site

Non-periodic lattice

lat = BuildLattice["Square", Bounded->{Region->Function[{r},Norm[r]<5]}]
sys =BuildSystem[lat,
HoppingRange->1,
Onsite->Function[{r},2],
Hopping->Function[{r1, r2},1]]

PlotLattice[lat]

PlotLattice[sys]

PlotLattice[sys,Labels->"Numbers", SiteTooltips->True, LinkTooltips->True]

To access the Hamiltonian matrix use sys[“H”]

sys["H"]

MatrixForm[sys["H"]]

sys["H"][[23,23]]
sys["H"][[23,24]]

Periodic lattice

lat = BuildLattice["Linear", Supercell->3]
sys =BuildSystem[lat,
HoppingRange->1,
Onsite->Function[{r},2],
Hopping->Function[{r1, r2},-1]]

PlotLattice[lat]

PlotLattice[sys, BoxedCells->False]

PlotLattice[sys,Labels->"Numbers", SiteTooltips->True, LinkTooltips->True]

sys[“H”] in this case gives the intra-unitcell Hamiltonian

MatrixForm[sys["H"]]

sys[“HBloch”] gives also the inter-unitcell Hamiltonians

sys["HBloch"]

MatrixForm/@sys["HBloch"][[All,2]]

{1}/.sys["HBloch"]
%[[1,3]]

Limitation: a unit cell can only have hoppings to nearest-neighbor unit cells (including diagonals). If longer range hopping is required, just enlarge the unitcell with Supercell.

Several orbitals per site

We just have to use a matrix Onsite and Hopping functions

lat = BuildLattice["Linear", Supercell->3]
sys =BuildSystem[lat,
HoppingRange->1,
Onsite->Function[{r},(({{2, 0}, {0, 2}}))],
Hopping->Function[{r1, r2},(({{-1, 0}, {0, -1}}))]]

PlotLattice[sys, SiteTooltips->True, LinkTooltips->True]

MatrixForm[sys["H"]]
MatrixForm[sys["Flat"]["H"]]

Example: Kane-Mele model (spin 1/2 per site)

lat = BuildLattice["Honeycomb", Tags->{"A","B"}]

With[{λ=0.1, t=3.16},KaneMele=Function[{r1,r2,tag1,tag2},If[Norm[r1-r2]≤1/(Sqrt[3]),t(({{1, 0}, {0, 1}})),I λ (({{1, 0}, {0, -1}}))Cos[3ArcTan@@(r1-r2)]If[tag1==="A",1,-1]]]];
sys =BuildSystem[lat,HoppingRange->1,Onsite->Function[{r},(({{0, 0}, {0, 0}}))],Hopping->KaneMele]

PlotLattice[sys, SiteTooltips->True, LinkTooltips->True, BoxCenterAtOrigin->True]

Example: BCS superconductor (Nambu representation H = (({{(Subscript[H, 0]), (SuperPlus[Δ])}, {Δ, -(Subscript[H, 0]^( * ))}})) )

lat = BuildLattice["Cubic"]
With[{t=1, Δ=0.1, μ=0.4},
sys =BuildSystem[lat,HoppingRange->1,
Onsite->Function[{r},(({{4t-μ, 0, 0, Δ}, {0, 4t-μ, -Δ, 0}, {0, -Δ, -(4t-μ), 0}, {Δ, 0, 0, -(4t-μ)}}))],
Hopping->Function[{r1, r2},(({{-t, 0, 0, 0}, {0, -t, 0, 0}, {0, 0, t, 0}, {0, 0, 0, t}}))],
Nambu->True]
]

PlotLattice[sys, SiteRadius->.3,LinkRadius->.08,SiteTooltips->True, LinkTooltips->True,BoxCenterAtOrigin->True]

Limitation: currently, all sites should have the same “SiteDimensions” (number of orbitals per site). If they are different, one should add “dummy” decoupled orbitals

Lesson 3 - Computing electronic structure

Non-periodic systems

Let’s build a spinless px+ipy 2D superconductor model on a disk. But this time we make sys a function of a parameter (the gap)

diskregion=Function[{r},Norm[r]≤25];
lat = BuildLattice["Square",  Bounded->{Region->diskregion}]
k[r1_,r2_]=-I ((r1-r2)/((Norm[r1-r2])^2));
With[{t=1.0, μ=1.0},
pwavehoppingΔ[Δ_]=Function[{r1,r2},(({{-t, 0}, {0, t}}))+Δ(({{0, k[r1,r2].{1,-I}}, {k[r1,r2].{1,I}, 0}}))];
sysΔ[Δ_] =BuildSystem[lat,HoppingRange->1,Onsite->Function[{r},(({{-μ, 0}, {0, μ}}))],Hopping->pwavehoppingΔ[Δ]]];

PlotLattice[sysΔ[.1]]

eigs0=Spectrum[sysΔ[0]];
eigs = Spectrum[sysΔ[0.1]];

ListPlot[{eigs0, eigs},PlotStyle->{Black, Red}, PlotDefaults, FrameLabel->{"eigenvalue","energy"}, BaseStyle->PointSize[.001]]
SmoothHistogram[{eigs0, eigs},.07,PlotStyle->{Black, Red}, PlotDefaults, FrameLabel->{"energy","DOS"}]

eigs0=Spectrum[sysΔ[0], Levels->800];
eigs = Spectrum[sysΔ[0.5],Levels->800];

ListPlot[{eigs0, eigs}, PlotStyle->{Black, Red},PlotDefaults,FrameLabel->{"eigenvalue","energy"}]

Eigenstates can also be computed, preserving their site/orbital structure. The orbital structure can also be processed with the OrbitalFunction option.

Eigenstates[sysΔ[0.5],Levels->2]
Dimensions[%]

Eigenstates[sysΔ[0.5],Levels->2, OrbitalFunction->Function[{ψ},((Norm[ψ])^2)]];
Dimensions[%]

We may visualize eigenstates on the lattice

PlotEigenstates[sysΔ[0.5],EnergyOrigin->0.2,Levels->4,PrintEigenvalues->True,MaskRegion->diskregion]

PlotEigenstates[sysΔ[0.5],EnergyOrigin->1.0,Levels->4,PrintEigenvalues->True, MaskRegion->diskregion]

We may also choose to plot a given function of the orbitals on each sites. For example, only the electron-like part of the wavefunction

PlotEigenstates[sysΔ[0.5],EnergyOrigin->1.0,Levels->4,PrintEigenvalues->True,MaskRegion->diskregion, OrbitalFunction->Function[{ψ},((Abs[ψ[[1]]])^2)]]

A rather slow “BubblePlot” method is also available

PlotEigenstates[sysΔ[0.5],EnergyOrigin->0.2,Levels->4,PrintEigenvalues->True,BubblePlot->True]

PlotEigenstates[sysΔ[0.5],EnergyOrigin->1.0,Levels->4,PrintEigenvalues->True,BubblePlot->True]

Kernel Polynomial Method for the DOS

ListLinePlot[DOSKPM[sysΔ[0.1], RandomSamples->4,Rescale->5.0, Order->300]]

The KPM accuracy grows with system size. It can compute spectral densities for systems way beyond the reach of exact diagonalization

diskregion=Function[{r},Norm[r]≤100];
lat = BuildLattice["Square",  Bounded->{Region->diskregion}]
k[r1_,r2_]=-I ((r1-r2)/((Norm[r1-r2])^2));
With[{t=1.0, μ=1.0},
pwavehoppingΔ[Δ_]=Function[{r1,r2},(({{-t, 0}, {0, t}}))+Δ(({{0, k[r1,r2].{1,-I}}, {k[r1,r2].{1,I}, 0}}))];
sysΔ[Δ_] =BuildSystem[lat,HoppingRange->1,Onsite->Function[{r},(({{-μ, 0}, {0, μ}}))],Hopping->pwavehoppingΔ[Δ]]];

ListLinePlot[DOSKPM[sysΔ[0.1], RandomSamples->4,Rescale->5.1, Order->300]]

Periodic systems

Quasi-1D system

A px+ipy superconducting nanoribbon

ribbonregion=Function[{r},0<=r[[1]]≤10];
lat = BuildLattice["Square", Bounded->{Axes->{1},Region->ribbonregion}]
k[r1_,r2_]=-I ((r1-r2)/((Norm[r1-r2])^2));
With[{t=1, μ=0.5, Δ=0.5},
pwavehopping=Function[{r1,r2},(({{-t, 0}, {0, t}}))+Δ(({{0, k[r1,r2].{1,-I}}, {k[r1,r2].{1,I}, 0}}))];
sys=BuildSystem[lat,HoppingRange->1,Onsite->Function[{r},(({{-μ, 0}, {0, μ}}))],Hopping->pwavehopping]];

PlotLattice[sys, BoxedCells->False]

We want to obtain the bandstructure of this periodic system using Bloch’s theorem: 
	Eigenstates of wavenumber k can be obtained by diagonalising H(k)=(Subscript[H, 0])+(∑_j)(e^(i k·(Subscript[A, j])))(Subscript[V, j]), where (Subscript[H, 0]) is the intra-cell and (Subscript[V, j]) the inter-cell Hamiltonians

To build the k Bloch harmonic of a periodic system sys, use sys[(Subscript[ϕ, 1]),(Subscript[ϕ, 2]),...] or sys[{(Subscript[ϕ, 1]),(Subscript[ϕ, 2]),...}], where (Subscript[ϕ, j])=k·(Subscript[A, j]). The full Brillouin zone is spanned by (Subscript[ϕ, j])∈[0,2π)

sys
sys[ϕ]

MatrixForm[sys["H"]]
MatrixForm[sys[ϕ]["H"]]

This is manual way to obtain the bandstructure of the nanoribbon

bands=Table[Spectrum[sys[ϕ]],{ϕ,0,2π,2π/120}];
ListPlot[Transpose[bands], DataRange->{0,2π}]

Given how often one needs to do this, we have two custom commands

bands=Bands[sys[ϕ],{ϕ,0,2π}, SamplePoints->120];
PlotBands[bands]

An more sophisticated UntangledBands algorithm can be used to connect adjacent eigenvalues correctly (using eigenstates)

bands=UntangledBands[sys[ϕ],{ϕ,0,2π},SamplePoints->120];
PlotBands[bands]

PlotEigenstates[sys[1.01π],Levels->2]
PlotEigenstates[sys[1.01π],Levels->2, BubblePlot->True, Overlay->True]
ListLinePlot[
Eigenstates[sys[1.01π],Levels->2, OrbitalFunction->Function[{ψ},((Norm[ψ])^2)]],
FrameLabel->{"site","|ψ((|)^2)"},PlotDefaults]

2D system

2D chiral superconductor

ribbonregion=Function[{r},0<=r[[1]]≤10];
lat = BuildLattice["Square"]
k[r1_,r2_]=-I ((r1-r2)/((Norm[r1-r2])^2));
With[{t=1, μ=0.5, Δ=0.5},
pwavehopping=Function[{r1,r2},(({{-t, 0}, {0, t}}))+Δ(({{0, k[r1,r2].{1,-I}}, {k[r1,r2].{1,I}, 0}}))];
sys=BuildSystem[lat,HoppingRange->1,Onsite->Function[{r},(({{-μ, 0}, {0, μ}}))],Hopping->pwavehopping]];

bands=Bands[sys[ϕ1, ϕ2],{ϕ1,0,2π},{ϕ2,0,2π}];
PlotBands[bands, AxesLabel->{"(Subscript[ϕ, 1])","(Subscript[ϕ, 2])","energy"}]

Graphene

sys=BuildSystem[BuildLattice["Honeycomb"],HoppingRange->1/(Sqrt[3])]
PlotBands[Bands[sys[{kx,ky}.(({{Cos[π/3], -Cos[π/3]}, {Sin[π/3], Sin[π/3]}}))],{kx,-((4π)/3),((4π)/3)},{ky,-((2π)/(Sqrt[3])),((2π)/(Sqrt[3]))}],AxesLabel->{"(Subscript[k, x])","(Subscript[k, y])","energy"}]

The Brillouin zone extends periodically beyond the first

PlotBands[Bands[sys[{kx,ky}.(({{Cos[π/3], -Cos[π/3]}, {Sin[π/3], Sin[π/3]}}))],{kx,-3((4π)/3),3((4π)/3)},{ky,-3((2π)/(Sqrt[3])),3((2π)/(Sqrt[3]))}],AxesLabel->{"(Subscript[k, x])","(Subscript[k, y])","energy"},
MeshFunctions->{Function[{x,y,z},y]}, Mesh->1]

A 1D cut of the above at (Subscript[k, y])=0

PlotBands[Bands[sys[{kx,0}.(({{Cos[π/3], -Cos[π/3]}, {Sin[π/3], Sin[π/3]}}))],{kx,0,4π}], FrameLabel->{"(Subscript[k, x])","energy"}]
PlotBands[UntangledBands[sys[{kx,0}.(({{Cos[π/3], -Cos[π/3]}, {Sin[π/3], Sin[π/3]}}))],{kx,0,4π}],FrameLabel->{"(Subscript[k, x])","energy"}]

If we want to “wrap” one dimension only, use the keyword None for the dimension that should remain open

PlotLattice[sys, UnitCells->3]

PlotLattice[sys[ϕ, None],UnitCells->3,BoxedCells->False]
PlotLattice[sys[None,ϕ],UnitCells->3,BoxedCells->False]

Missing Bloch phases are interpreted as trailing None’s

sys[ϕ]==sys[ϕ,None]

Lesson 4 - Building a ScatteringSystem

A transport system consists of N semi-infinite periodic leads, attached to a central system.

lead1=BuildLattice["Honeycomb",Tags->{"1"},Bounded->{Axes->{1},Region->Function[{r},Abs[r[[1]]]<2]}]
lead2=BuildLattice["Honeycomb",Tags->{"2"},Bounded->{Axes->{2},Region->Function[{r},Abs[r[[1]]]<3]}]

central = BuildLattice["Triangular",Tags->{"C"},Bounded->{Axes->All,Region->Function[{r},Norm[r]<9]}]

{PlotLattice[lead1],PlotLattice[lead2],PlotLattice[central]}

To attach the leads to the central region, we need to reshape their unit cell to the shape of the central region (its complementary).

lead1=BuildLattice[lead1, Reshape->{Region->Function[{r},Norm[r]>9]}]
lead2=BuildLattice[lead2, Reshape->{Region->Function[{r},Norm[r]>9]}]

PlotLattice[lead1]

PlotLattice[{lead1,lead2,central}, UnitCells->10, BoxedCells->False, SiteColors->{{Red},{Green},{Blue}}, Labels->"Tags"]

The scattering system is built just like a regular system, but using BuildScatteringSystem[central, {lead1,lead2,...}, opts...]
A useful option is CouplingFactor, which multiplies the hoppings between leads and central region by a constact factor.

scatsys=BuildScatteringSystem[central,{lead1,lead2},HoppingRange->1,
Hopping->Function[{r1,r2,tag1,tag2},If[tag1==tag2,If[tag1=="C",1,If[Norm[r1-r2]>1/(Sqrt[3]),0,1]],1]],
CouplingFactor->0.5]

PlotLattice[scatsys,UnitCells->10, BoxedCells->False, SiteColors->{{Darker[Red]},{Darker[Green]},{Darker[Blue]}}, Labels->"Numbers"]

scatsys["Central"]
scatsys["Leads"]
scatsys["Couplings"]
scatsys["Nambu"]

Lesson 5 - Computing transport

Normal systems

Let us rebuild our sample scattering system with a magnetic flux in the central region (Peierls phase in the hoppings)

lead1=BuildLattice["Honeycomb",Tags->{"1"},Bounded->{Axes->{1},Region->Function[{r},Abs[r[[1]]]<2]},Reshape->{Region->Function[{r},Norm[r]>9]}]
lead2=BuildLattice["Honeycomb",Tags->{"2"},Bounded->{Axes->{2},Region->Function[{r},Abs[r[[1]]]<3]},Reshape->{Region->Function[{r},Norm[r]>9]}]
central = BuildLattice["Triangular",Tags->{"C"},Bounded->{Axes->All,Region->Function[{r},Norm[r]<9]}]

With[{A=Function[{r},{0.,r[[1]]}]},
scatsysB[B_]=BuildScatteringSystem[central,{lead1,lead2},HoppingRange->1,
Hopping->Function[{r1,r2,tag1,tag2},If[tag1==tag2,If[tag1=="C",(E^(I B A[((r1+r2)/2)].(r1-r2))),If[Norm[r1-r2]>1/(Sqrt[3]),0,1]],1]]]]

PlotLattice[scatsysB[0]]

We may compute the bandstructure of the leads

PlotBands[Bands[scatsysB[0]["Leads"][[1]][ϕ],{ϕ,0,2π}]]
PlotBands[Bands[scatsysB[0]["Leads"][[2]][ϕ],{ϕ,0,2π}]]

We compute the conductance as a function of Fermi energy at zero flux.

Table[{ϵF,Conductance[ϵF, scatsysB[0]]},{ϵF, 0.1,0.5,0.001}];
ListLinePlot[%, PlotDefaults]

Sampling can also be done using Sample. We do a finite flux.

Sample[Conductance[ϵF, scatsysB[.8]],{ϵF,0.1,0.5}, SamplePoints->300];
ListLinePlot[%, PlotDefaults]

We can compute the full scattering matrix, S(ω)=(({{(Subscript[t, 11]), (Subscript[t, 12])}, {(Subscript[t, 21]), (Subscript[t, 22])}})). Each (Subscript[t, ij]) submatrix has a dimension equal to the number of propagating modes at a given energy ω in lead i and j.

Smatrix[0.4,scatsysB[.8]]
MatrixForm@Map[MatrixForm,%,{2}]

MatrixForm[Smatrix[2,2,0.4,scatsysB[.8]]]

We can compute scattering states

ScatteringStates[0.4, scatsysB[0]]

We can plot scattering states incoming from each of the leads at a given energy

With[{B=0, ϵF=0.4},
PlotScatteringStates[scatsysB[B],
ScatteringStates[ϵF, scatsysB[B],LeadCells->25],
Lead->1,BubblePlot->True]]

With[{B=0, ϵF=0.4},
PlotScatteringStates[scatsysB[B],
ScatteringStates[ϵF, scatsysB[B],LeadCells->25],
Lead->2,BubblePlot->True]]

With[{B=0.8, ϵF=0.4},
PlotScatteringStates[scatsysB[B],
ScatteringStates[ϵF, scatsysB[B],LeadCells->25],
Lead->2,BubblePlot->True]]

Superconducting systems

NS systems

We can build a normal lead coupled to a (grounded) superconducting island

central = BuildLattice["Square",Tags->"S",Bounded->{Axes->All,Region->Function[{r},Abs[r[[1]]]<6 &&Abs[r[[2]]]<6]}]
lead=BuildLattice["Square",Tags->"N",Supercell->-1,Bounded->{Axes->{2},Region->Function[{r},Abs[r[[2]]]<3]},
Reshape->{Region->Function[{r},r[[1]]≤-6]}]

With[{μ=0.2, Δ=0.5, t=0.1},
scatsys=BuildScatteringSystem[central,{lead},
Onsite->Function[{r,tag},If[tag=="N",(({{-μ, 0}, {0, μ}})),(({{-μ, Δ}, {Δ, μ}}))]],
Hopping->Function[{r1, r2},(({{-t, 0}, {0, t}}))],
Nambu->True]
]

PlotLattice[scatsys, UnitCells->5, SiteTooltips->True]

PlotBands[UntangledBands[scatsys["Leads"][[1]][ϕ],{ϕ,0,2π}]]

While the system has zero normal conductance, there is a finite BTK conductance. The Nambu->True option indicates that BTK should be used

Sample[Conductance[1,ϵF,scatsys],{ϵF,-0.5,0.5}, SamplePoints->300];
ListLinePlot[%,PlotRange->{0,0.5},PlotDefaults]

states=ScatteringStates[0.045, scatsys, LeadCells->20]
PlotScatteringStates[scatsys,states,PlotFunction->Function[{ψ},((Abs[ψ[[1]]])^2)],PlotLegends->False]
PlotScatteringStates[scatsys,states,PlotFunction->Norm,BubblePlot->True]

Lesson 6 - Advanced topics [Work in progress]

Supercurrents

BuildLattice["Linear",Bounded->{Region->(Abs[#[[1]]]<100&)}]

With[{Δ=0.2, t=1, μ=0.2},
sysϕ[ϕ_, B_]=BuildSystem[BuildLattice["Linear",Bounded->{Region->(Abs[#[[1]]]<100&)}],Onsite->((({{-μ, B, 0, Δ (E^(-I ϕ Sign[#[[1]]]/2 ))}, {B, -μ, -Δ (E^(-I ϕ Sign[#[[1]]]/2 )), 0}, {0, -Δ (E^(I ϕ Sign[#[[1]]]/2 )), μ, -B}, {Δ (E^(I ϕ Sign[#[[1]]]/2 )), 0, -B, μ}}))&), Hopping->(-t(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}}))&), Nambu->True]]

PlotBands[UntangledBands[sysϕ[ϕ, 0.03],{ϕ,0,2π}, Levels->22]]

JosephsonCurrent[sysϕ[ϕ,0.03],ϕ, Levels->22, SamplePoints->100];
ListLinePlot[%,PlotDefaults, FrameLabel->{"ϕ", "(Subscript[I, J])[e/ℏ]"}]
Max[%%[[All,2]]]

Elastic Relaxation: SpringMesh, ElasticEnergy, ElasticForces, RelaxLattice, GradientMesh, RelaxField

Band topology: MomentumCycles, BerryCurvature, BerryMesh, ChernNumber

Linear Response: LinearResponseIm, OpticalConductivity, LocalDensityOfStates