Condensed Matter Physics

Part 1. Atomic System

Crystal Geometry

Crystal Formation

Atomic Orders

Atomic order: how atoms are arranged in the matter.

Gas (no order) - atoms are well separated, no regular arrangement.

Liquid (short-ranged order) - atoms are close together, no regular arrangement.

Solid (long-ranged order) - atoms are tightly packed in a regular pattern (in particular, atoms form a regular lattice in a crystal).

The atomic order is characterized by the radial density distribution ρ(r) - density ρ of atoms as a function of the distance r from a reference atom (in unit of average inter-atomic spacing a).

Row[{Plot[Exp[-(1/r^12-1/r^6)/5],{r,0,4},Filling->Axis,PlotStyle->Black,FrameLabel->{"\!\(r/a\)","\!\(ρ(r)/ρ\_∞\)"},GridLines->{None,{1}},PlotRange->{0,3},PlotRangePadding->Scaled[.02],PlotLabel->"gas",AspectRatio->1,ImageSize->140],ListLinePlot[{{0.025815813728932546`,0.00008577913851359753`},{0.07475915115870785`,0.00008577913851359753`},{0.12370248858848315`,0.00008577913851359753`},{0.17264582601825834`,0.00008577913851359753`},{0.22158916344803364`,0.00008577913851359753`},{0.27053250087780895`,0.00008577913851359753`},{0.31947583830758425`,0.00008577913851359753`},{0.36841917573735944`,0.00008577913851359753`},{0.41736251316713474`,0.00008577913851359753`},{0.46630585059691004`,0.00008577913851359753`},{0.5152491880266854`,0.00008577913851359753`},{0.5641925254564606`,0.00008577913851359753`},{0.6131358628862358`,0.00008577913851359753`},{0.662079200316011`,0.00008577913851359753`},{0.7110225377457864`,0.00008577913851359753`},{0.7599658751755616`,0.00008577913851359753`},{0.808909212605337`,0.00008577913851359753`},{0.8578525500351122`,0.00008577913851359753`},{0.9067958874648874`,0.006653281921836562`},{0.9407225190923454`,0.05460213326119101`},{0.9506542028240366`,0.14342413353643835`},{0.9601886192064606`,0.23186389819995767`},{0.9646380135182582`,0.3256231084600225`},{0.969087407830056`,0.399292652027027`},{0.9745255564333644`,0.48794583157376126`},{0.9779861964536516`,0.5817719330658782`},{0.9824355907654492`,0.6823763460726351`},{0.9735368021418538`,0.6288633604307434`},{0.986884985077247`,0.8164932663376265`},{0.9846602879213482`,0.749267578125`},{0.9923231336805554`,0.9439953869885511`},{0.9824355907654492`,0.8803743929476351`},{0.9957837737008426`,1.0833892822265625`},{0.9913343793890448`,1.0168325063344597`},{1.0002331680126402`,1.2165028340107686`},{0.9980084708567414`,1.1479393211570947`},{1.004682562324438`,1.3663391938080658`},{1.001716299449906`,1.2924243823902026`},{1.0091319566362358`,1.5158410974451013`},{1.0061656937617038`,1.4422607421875`},{1.0145701052395442`,1.6699882255302176`},{1.0106150880735016`,1.5903133357967343`},{1.0191430938377808`,1.8205306594436232`},{1.0135813509480336`,1.7490685265343466`},{1.022480139571629`,1.9666880014780406`},{1.0180307452598312`,1.893553587767455`},{1.0269295338834268`,2.1158554489548145`},{1.025446402446161`,2.038038649000563`},{1.0313789281952246`,2.2526480185019`},{1.0291542310393256`,2.188766891891892`},{1.0364639502658506`,2.383324818261342`},{1.0291542310393256`,2.3252250052787162`},{1.04027771681882`,2.4994766647751265`},{1.0313789281952246`,2.4429535736908785`},{1.0457158654221284`,2.6106275953688063`},{1.0506596368796814`,2.7047212617891327`},{1.0536258997542134`,2.7725935652449323`},{1.057333728347378`,2.843409082911036`},{1.0749829924508425`,2.9304568728885134`},{1.0907041856858615`,2.986110377956081`},{1.1159174201193818`,2.883097880595439`},{1.1255577744616103`,2.799706811303491`},{1.1301554819171347`,2.7265723975929053`},{1.1381643916783706`,2.651654217694257`},{1.1470631803019662`,2.5516518257759713`},{1.1559619689255616`,2.4449603106524496`},{1.1648607575491572`,2.3345897777660474`},{1.1693101518609548`,2.2545878642314188`},{1.1776527661955756`,2.169368434596706`},{1.1826583347963482`,2.092978647592905`},{1.1915571234199436`,2.000133617504223`},{1.2004559120435392`,1.8978855913670367`},{1.2057951852176965`,1.8238083298141892`},{1.2138040949789324`,1.7438446398407335`},{1.2191433681530897`,1.673436840160473`},{1.2271522779143256`,1.59443829510663`},{1.2360510665379212`,1.503131763355152`},{1.2456914208801497`,1.4284365542300113`},{1.2538486437851122`,1.3603858741554054`},{1.2620058666900746`,1.2910865577491555`},{1.27090465531367`,1.2206277933206644`},{1.2798034439372656`,1.1559662690033785`},{1.2919863093148072`,1.08792833011583`},{1.3072413755266852`,1.0128190324113175`},{1.3205895584620784`,0.9490493911880629`},{1.3361624385533706`,0.8867290599926099`},{1.3556903358107053`,0.8212128143768771`},{1.3791732502340823`,0.7626458245354728`},{1.4095774446980334`,0.6984302417652026`},{1.4451725991924156`,0.6395659575591215`},{1.489666542310393`,0.5935934380758598`},{1.5386098797401684`,0.5736476888820641`},{1.5875532171699436`,0.583620563478962`},{1.6364965545997188`,0.617917522458538`},{1.6809904977176964`,0.6737548094969972`},{1.7165856522120784`,0.7320669755972489`},{1.7455067152387636`,0.7885104342623874`},{1.7722030811095504`,0.8469287769214526`},{1.7988994469803368`,0.9062390026745493`},{1.8255958128511232`,0.9673329946157097`},{1.8545168758778088`,1.03059298835666`},{1.8856626360603927`,1.0902217437861967`},{1.919033093398876`,1.1482737773173564`},{1.959077642205056`,1.205733345650338`},{2.0057962824789324`,1.2603165910050675`},{2.0547396199087076`,1.2938838274531632`},{2.1036829573384828`,1.2919379007025493`},{2.152626294768258`,1.253019365690264`},{2.1971202378862356`,1.19580371387012`},{2.2349400895365163`,1.1345610747466215`},{2.2683105468749996`,1.0772257329874515`},{2.3016810042134828`,1.0171669624947213`},{2.3395008558637636`,0.9573736333990239`},{2.383994798981741`,0.9039687547988331`},{2.4329381364115163`,0.8609151254414926`},{2.481881473841292`,0.8331856692452395`},{2.530824811271067`,0.8241857580236487`},{2.5797681487008424`,0.8339153917767197`},{2.6287114861306176`,0.8628610521921072`},{2.6776548235603927`,0.9068876449247543`},{2.726598160990168`,0.9567520179092441`},{2.7755414984199436`,1.0058866683622543`},{2.8244848358497188`,1.0486970568757679`},{2.873428173279494`,1.0849399426059585`},{2.922371510709269`,1.111696435426904`},{2.9713148481390443`,1.118507179054054`},{3.02025818556882`,1.1102369903639433`},{3.069201522998595`,1.0912642045454546`},{3.1181448604283704`,1.063777989193028`},{3.1670881978581455`,1.0321566794955466`},{3.2160315352879207`,1.0000488881104115`},{3.2649748727176964`,0.9676978558814495`},{3.3139182101474716`,0.9353468236524876`},{3.3628615475772468`,0.9112659801136362`},{3.411804885007022`,0.9000769012976044`},{3.460748222436797`,0.9076173674562344`},{3.5096915598665723`,0.9299955250882985`},{3.558634897296348`,0.9623465573172605`},{3.607578234726123`,1.0000488881104115`},{3.6565215721558983`,1.0382377005912162`},{3.7054649095856735`,1.0735076229460994`},{3.7544082470154487`,1.0881020735757065`},{3.803351584445224`,1.0812913299485567`},{3.852294921874999`,1.0557510413467441`},{3.901238259304775`,1.020481118991861`},{3.9501815967345504`,0.980589620604269`},{3.988001448384831`,0.9512790989231417`}},Filling->Axis,PlotStyle->Black,FrameLabel->{"\!\(r/a\)","\!\(ρ(r)/ρ\_∞\)"},GridLines->{None,{1}},PlotRange->{0,3},PlotRangePadding->Scaled[.02],PlotLabel->"liquid",AspectRatio->1,ImageSize->140],
Plot[Evaluate[12Total[#2 Exp[-((r-1.122#1)/0.03)^2/2]&@@@SortBy[Tally[Cases[Norm/@Tuples[Range[-4,4],3],x_/;0<x<4]],First[N[#]]&]]/(4π r^2)],{r,0,4},Filling->Axis,PlotStyle->Black,FrameLabel->{"\!\(r/a\)","\!\(ρ(r)/ρ\_∞\)"},GridLines->{None,{1}},PlotRange->{0,3},PlotRangePadding->Scaled[.02],PlotLabel->"solid",AspectRatio->1,ImageSize->140]}]

Many macroscopic behaviors can be explained

({{ , shape, volume, fluidity, compressibility}, {gas, ×, ×, high, high}, {liquid, ×, ✓, high, low}, {solid, ✓, ✓, low, low}})

Atomic Interaction

Order always arise from interaction.

Atoms (or molecules) interact in many different ways (i.e. different chemical bounds):

Ionic bound.

Block[{atom},
atom[col_,txt_,es_]:=Block[{r1=0.45,r2=0.7,r3=0.08},{{EdgeForm[Black],Lighter[col,0.8],Disk[{0,0},r1]},Text[Style[StringTemplate["\!\(`1`\)"][txt],SingleLetterItalics->False]],
If[Length[es]==0,{},{Circle[{0,0},r2]}],{EdgeForm[Black],Yellow,Disk[r2 ReIm@#,r3]&/@Exp[I π/4 es]}}];
Graphics[{Arrowheads[{{0.01,1,Graphics[Line[Offset/@{{-4,-3},{0,0},{-4,3}}]]}}],Translate[{Translate[atom[Gray,"Na",{0}],{-1,0}],Text["+"],{Red,Arrow[BSplineCurve[{{-0.2,0.2},{0,0.4},{0.2,0.2}}]]},Translate[atom[Gray,"Cl",{-0.25,0.25,1.75,2.25,4,5.75,6.25}],{1,0}]},{-2,0}],Text["→"],
Translate[{Translate[atom[Red,"Na\^+",{}],{-1,0}],{AbsoluteThickness[1],Translate[Line[{0,#}&/@{-0.1,0.1}],{#,0}]&/@Range[-0.35,0.1,0.2]},Translate[atom[Blue,"Cl\^-",{-0.25,0.25,1.75,2.25,3.75,4.25,5.75,6.25}],{1,0}]},{1.8,0}]},PlotRange->{All,{-1,1}},ImageSize->{Automatic,60}]]

Covalent bound.

Block[{atom},
atom[col_,txt_,es_]:=Block[{r1=0.45,r2=0.7,r3=0.08},{{EdgeForm[Black],Lighter[col,0.8],Disk[{0,0},r1]},Text[Style[StringTemplate["\!\(`1`\)"][txt],SingleLetterItalics->False]],
{Circle[{0,0},r2]},{EdgeForm[Black],Yellow,Disk[r2 ReIm@#,r3]&/@Exp[I π/4 es]}}];
Graphics[{Arrowheads[{{0.01,1,Graphics[Line[Offset/@{{-4,-3},{0,0},{-4,3}}]]}}],Translate[{Translate[atom[Gray,"H",{0}],{-1,0}],Text["+"],{Red,Rotate[Arrow[BSplineCurve[{{-0.2,0.2},{0,0.4},{0.2,0.2}}]],#,{0,0}]&/@{0,π}},Translate[atom[Gray,"H",{4}],{1,0}]},{-2,0}],Text["→"],
Translate[{Translate[atom[Gray,"H",{}],{-0.7,0}],Translate[atom[Gray,"H",{3.75,4.25}],{0.7,0}]},{1.7,0}]},PlotRange->{All,{-1,1}},ImageSize->{Automatic,60}]]

Metallic bound.

Block[{r1=0.45,r2=0.7,r3=0.08,atom},
atom[col_,txt_,es_]:={{EdgeForm[Black],Lighter[col,0.8],Disk[{0,0},r1]},Text[Style[StringTemplate["\!\(`1`\)"][txt],SingleLetterItalics->False]],
If[Length[es]==0,{},{Circle[{0,0},r2]}],{EdgeForm[Black],Yellow,Disk[r2 ReIm@#,r3]&/@Exp[I π/4 es]}};
Graphics[{Arrowheads[{{0.01,1,Graphics[Line[Offset/@{{-4,-3},{0,0},{-4,3}}]]}}],Translate[{Translate[atom[Gray,"Na",{2}],{-1,0}],Text["+"],Translate[atom[Gray,"Na",{-2}],{1,0}],Text["+",{2,0}],Text["…",{2,0},{-2,0}]},{-3,0}],Text["→"],
Block[{x1,x2,softmax},
{x1,x2}={-0.6,0.6};
softmax[ys_]:=(#.ys/Total[#])&@Exp[4ys];
Translate[{Translate[atom[Red,"Na\^+",{}],{x1,0}],Translate[atom[Red,"Na\^+",{}],{x2,0}],Text["…",{-1.5,0}],Text["…",{1.5,0}],{Translate[Rotate[Line[Table[{x,softmax[Re[#]-Im[#]^2&@Sqrt[r2^2-(x-{x1,x2})^2]]},{x,x1,x2,(x2-x1)/100}]],#,{0,0}]&/@{0,π},{#,0}]&/@{-2x2,0,2x2},
{EdgeForm[Black],Yellow,Disk[#,r3]&/@{{x1,r2},{x2,-r2}}}}},{2,0}]]},PlotRange->{All,{-1,1}},ImageSize->{Automatic,60}]]

Hydrogen bound.

Block[{r1=0.45,r2=0.7,r3=0.08,atom,H2O},
atom[col_,txt_,es_]:={{EdgeForm[Black],Lighter[col,0.8],Disk[{0,0},r1]},Text[Style[StringTemplate["\!\(`1`\)"][txt],SingleLetterItalics->False]],
{Circle[{0,0},r2]},{EdgeForm[Black],Yellow,Disk[r2 ReIm@#,r3]&/@Exp[I π/4 es]}};
H2O={Translate[atom[Red,"H",{}],2r2{1,-1}/Sqrt[2]],Translate[atom[Red,"H",{}],2r2{-1,-1}/Sqrt[2]],atom[Blue,"O",{0.75,1.25,2.75,3.25,4.75,5.25,6.75,7.25}]};
Graphics[{Translate[H2O,{-1.5,1.5}],Translate[H2O,{1.5,-1.5}],{Translate[{AbsoluteThickness[1],Rotate[Translate[Line[{0,#}&/@{-0.1,0.1}],{#,0}]&/@Range[-0.4,0.4,0.2],-π/4,{0,0}],Text["δ+",{0.,0.7}],Text["δ-",{0.,-0.7}]},{1,-1}/2]}},PlotRange->{All,{-3.5,2.5}},ImageSize->{Automatic,150}]]

Van der Waals force (Molecular bound).

Block[{mol,efield},
mol={First@ContourPlot[y,{x,-1,1},{y,-1,1},ColorFunction->(Lighter[ColorData["BlueRedSplit"][#]]&),ContourStyle->None,RegionFunction->(#1^2+#2^2<=1&),BoundaryStyle->Black],{Red,Text[Style["+",Bold],#]&/@ReIm[1.25Exp[I π/2(1+0.3 Range[-1,1])]]},{Blue,Text[Style["-",Bold],#]&/@ReIm[1.25Exp[-I π/2(1+0.3 Range[-1,1])]]}};
efield={First@ContourPlot[Re[1/(x+I y)],{x,-6,6},{y,-6,6},ContourLabels->None,ContourStyle->AbsoluteThickness[1],Contours->7,ContourShading->None],{AbsoluteThickness[1],Translate[Rotate[Line[Offset/@{{-4,3},{0,0},{-4,-3}}],#1,{0,0}],#2]&@@@{{π/2,{0,2}},{π/2,{0,-2}},{-π/2,{-1.96,0}},{-π/2,{1.96,0}},{-π/2,{2.94,0}},{0,{3,2.94}},{π,{3,-2.94}}}}};
Graphics[{Translate[efield,{-3,0}],Translate[mol,{-3,0}],Translate[Rotate[mol,π,{0,0}],{3,0}]},PlotRange->{{-5.2,5.2},{-3.2,3.2}},ImageSize->180]]

All types of atomic forces are the residual of the electrostatic force of electrons and nuclei.

They are generally repulsive at short range, and attractive at longer range. The pair-wise interaction can be described by the inter-atomic potential energy V(r) (where r is the inter-atom distance). One famous phenomenological model (the 6-12 potential) is

V(r)=(Subscript[V, 0])((1/(((r/a))^12))-(2/(((r/a))^6))).

Plot[1/r^12-2/r^6,{r,0,3},PlotStyle->Black,PlotRange->{-1.5,2},FrameLabel->{"\!\(r/a\)","\!\(V(r)/V\_0\)"},GridLines->{{1},None},Prolog->{{{Lighter[Red,0.8],Rectangle[{0,-2},{1,2}]},{Red,Text["repul.",{0.4,0.5}]}},{{Lighter[Blue,0.8],Rectangle[{1,-2},{3,2}]},{Blue,Text["attr.",{1.5,0.5}]}}},ImageSize->200]

They can be isotropic (ionic bounds, metallic bounds, Van der Waals) or directional (most covalent bounds, hydrogen bounds).

Statistical Mechanics of Atoms

Let (Subscript[r, i]) be the position of the ith atom.

The total potential energy of all atoms is

V((Subscript[r, 1]),(Subscript[r, 2]),…)=(∑_(i≠j))V((Subscript[r, i])-(Subscript[r, j])),

(summing the potential energy over all pairs of atoms).

For isotropic interaction, V(r) only depends on r=|r|, and is often denoted as V(r).

For directional interaction, V(r) also depends on the direction of r. [But we generally assume the inversion symmetry V(r)=V(-r), i.e. the action and reaction are always equal in magnitude and opposite in direction (Newton’s third law).]

By Boltzmann distribution, the probability (density) to observe the atoms at positions {(Subscript[r, i])} is given by

p((Subscript[r, 1]),(Subscript[r, 2]),…)=(1/Z)exp(-((V((Subscript[r, 1]),(Subscript[r, 2]),…))/((Subscript[k, B])T))),
Z=∫d(Subscript[r, 1])d(Subscript[r, 2])… exp(-((V((Subscript[r, 1]),(Subscript[r, 2]),…))/((Subscript[k, B])T))).

(Subscript[k, B]) - Boltzmann constant, T - temperature.

High-temperature limit (T→∞): p((Subscript[r, 1]),(Subscript[r, 2]),…)=1/Z. The atoms can appear anywhere with equal probability independently. => Gas phase.

Low-temperature limit (T->0): the probability is dominated by the configuration of the lowest energy. The atoms will arrange themselves to minimize the total potential energy

(min_({(Subscript[r, i])})) V((Subscript[r, 1]),(Subscript[r, 2]),…).

This leads to crystal formation. => Solid phase.

Minimal Energy Configurations

Isotropic Repulsive Interaction

V(r)=((Subscript[V, 0])/(((r/a))^12)).

Plot[1/(0.1+r^12),{r,0,3},PlotStyle->Black,PlotRange->{-1.5,2},FrameLabel->{"\!\(r/a\)","\!\(V(r)/V\_0\)"},ImageSize->200]

Block[{L=20,V,U,rs0,plt,U0,sol},
V=Function[{Typed[r,"Real64"]},(1)/(0.1+r^12)];
U[rs_List]:=Total[V[Norm[#1-#2]]&@@@Subsets[rs,{2}]];
rs0=Flatten[Table[{x,y}+RandomReal[{-0.2,0.2},2],{x,-2,2},{y,-2,2}],1];
plt[rs_List]:=Graphics[Point[rs],PlotRange->L/2,Frame->True,FrameTicks->None,ImageSize->140];
{U0,sol}=Monitor[FindMinimum[U[rs],{{rs,rs0}},Method->"QuasiNewton",AccuracyGoal->3,MaxIterations->1000],plt[rs]];plt[rs/.sol]]

Isotropic Attractive Interaction

V(r)=-((2(Subscript[V, 0]))/(((r/a))^6)).

Plot[-2/(0.5+r^6),{r,0,3},PlotStyle->Black,PlotRange->{-1.5,2},FrameLabel->{"\!\(r/a\)","\!\(V(r)/V\_0\)"},ImageSize->200]

Block[{L=20,V,U,rs0,plt,U0,sol},
V=Function[{Typed[r,"Real64"]},-2/(0.5+r^6)];
U[rs_List]:=Total[V[Norm[#1-#2]]&@@@Subsets[rs,{2}]];
rs0=Flatten[Table[{x,y}+RandomReal[{-0.2,0.2},2],{x,-2,2},{y,-2,2}],1];
plt[rs_List]:=Graphics[Point[rs],PlotRange->L/2,Frame->True,FrameTicks->None,ImageSize->140];
{U0,sol}=Monitor[FindMinimum[U[rs],{{rs,rs0}},Method->"QuasiNewton",AccuracyGoal->3,MaxIterations->1000],plt[rs]];plt[rs/.sol]]

Isotropic Interaction with Both Attraction and Repulsion

V(r)=(Subscript[V, 0])((1/(((r/a))^12))-(2/(((r/a))^6))).

Plot[(1-2r^6)/(0.1+r^12),{r,0,3},PlotStyle->Black,PlotRange->{-1.5,2},FrameLabel->{"\!\(r/a\)","\!\(V(r)/V\_0\)"},ImageSize->200]

Block[{L=20,V,U,rs0,plt,U0,sol},
V=Function[{Typed[r,"Real64"]},(1-2r^6)/(0.1+r^12)];
U[rs_List]:=Total[V[Norm[#1-#2]]&@@@Subsets[rs,{2}]];
rs0=Flatten[Table[{x,y}+RandomReal[{-0.2,0.2},2],{x,-2,2},{y,-2,2}],1];
plt[rs_List]:=Graphics[Point[rs],PlotRange->L/2,Frame->True,FrameTicks->None,ImageSize->140];
{U0,sol}=Monitor[FindMinimum[U[rs],{{rs,rs0}},Method->"QuasiNewton",AccuracyGoal->3,MaxIterations->1000],plt[rs]];plt[rs/.sol]]

Need short-range repulsion + long-range attraction to crystalize.
=> The existence of crystals implies that the inter-atomic interaction is neither purely repulsive (unlike electric force) nor purely attractive (unlike gravity).

We can infer microscopic laws of physics by observing their macroscopic consequences.

The crystal aways have roughly the same volume.

Because the atomic spacing is determined by the equilibrium position r=a of the inter-atomic potential (the distance where interaction switch from repulsive to attractive). => Explains why a solid is hard to compress.

The crystal aways assume the same triangular lattice structure.

Because triangular lattice is the closest packed structure in 2D.

Spontaneous symmetry breaking: although the microscopic interaction is isotropic (continuous rotation symmetry, no special bound angle), the resulting low-energy configuration (the crystal lattice) spontaneously breaks the continuous rotation symmetry to a discrete (6-fold) rotation symmetry with an emergent  60° bound angle. => Explains why a solid is rigid (has a fixed shape).

This is an emergent phenomenon. No where in the microscopic theory can we see why 60° bound angle should be special. => Demonstrates that reductionism doesn’t always work.

Directional Interaction

V(r)=(Subscript[V, 0])((1/(((r/a))^12))-((1+cos 4θ)/(((r/a))^6))),

Row@{Plot3D[Block[{r=Norm[{x,y}],θ=ArcTan[x,y]},(1-(1+Cos[4θ])r^6)/(0.1+r^12)],{x,-2,2},{y,-2,2},PlotRange->{-1,1},PlotPoints->32,ClippingStyle->None,ColorFunction->"BlueRedSplit",MeshFunctions->{#3&},AxesLabel->(Style[#,FontFamily->"LaTeX"]&/@{"\!\(,x\)","\!\(,y\)","\!\( V\)"}),TicksStyle->Directive[FontFamily->"LaTeX"],BoxRatios->{1, 1, 1},ViewPoint->{0.8,-2,0.4},ImageSize->180],
ContourPlot[Block[{r=Norm[{x,y}],θ=ArcTan[x,y]},(1-(1+Cos[4θ])r^6)/(0.1+r^12)],{x,-2,2},{y,-2,2},Contours->64,ContourStyle->None,PlotLabel->"top view",PlotRange->{-1,1},ClippingStyle->None,FrameLabel->{"\!\(,x\)","\!\(,y\)"},ImageSize->140]}

Block[{L=20,V,U,rs0,plt,U0,sol},
V[{x_,y_}]:=Block[{r=Norm[{x,y}],θ=ArcTan[x,y]},(1-(1+Cos[4θ])r^6)/(0.1+r^12)];
U[rs_List]:=Total[V[#1-#2]&@@@Subsets[rs,{2}]];
rs0=Flatten[Table[{x,y}+RandomReal[{-0.3,0.3},2],{x,-2,2},{y,-2,2}],1];
plt[rs_List]:=Graphics[Point[rs],PlotRange->L/2,Frame->True,FrameTicks->None,ImageSize->140];
{U0,sol}=Monitor[FindMinimum[U[rs],{{rs,rs0}},Method->"QuasiNewton",AccuracyGoal->3,MaxIterations->1000],plt[rs]];plt[rs/.sol]]

The explicit symmetry breaking in the microscopic inter-atomic potential already breaks the continuous rotation symmetry to 4-fold discrete rotation symmetry.

The square lattice is compatible with the 4-fold rotation symmetry and gains energy from bounds in all four directions.

The triangular lattice is incompatible with the 4-fold rotation symmetry, and can not gain as much energy, thus it ceases to be the minimal energy configuration.

Gas-Solid Transition

N: number of atoms in the system,

V: volume of the system,

T: temperature of the system (in equilibrium with a thermal bath).

F=E-T S,

E: total energy of the system,

S: entropy of the system,

S=(Subscript[k, B])ln Ω,

Ω: number of microscopic configurations

Gas phase: every atom can reach every point in the space, number of configurations for a single atom ∝V, number of configurations for N atoms ∝(V^N); however, permutations of atoms does not leads to new configurations, thus a factor N!  must be divided,

(Subscript[Ω, gas])=((V^N)/(N!))
=>(Subscript[S, gas])=(Subscript[k, B])(N ln V-ln N!)
≃(Subscript[k, B])(N ln V-N ln N+N)

Solid phase: every atom oscillates typically within the V/N volume surrounding its equilibrium position on the lattice, total number of configurations ∝(((V/N))^N). Moreover, atoms are distinguishable by their lattice coordinates, hence there is no permutation redundancy.

(Subscript[Ω, solid])=((((V/N)))^N)
=>(Subscript[S, solid])=(Subscript[k, B])(N ln V-N ln N)

(Subscript[F, gas])-(Subscript[F, solid])
=((Subscript[E, gas])-(Subscript[E, solid]))-T((Subscript[S, gas])-(Subscript[S, solid]))
=(Subscript[E, bind])-N (Subscript[k, B])T

(Subscript[E, bind]) is the binding energy, i.e. the energy released when atoms binds to solid. It defines a temperature scale

(SubStar[T])=((Subscript[E, bind])/(N (Subscript[k, B]))).

When T<(SubStar[T]), (Subscript[F, gas])>(Subscript[F, solid]), the system is in the solid phase,

When T>(SubStar[T]), (Subscript[F, gas])<(Subscript[F, solid]), the system is in the gas phase,

Therefore (SubStar[T]) is the gas-solid transition temperature.

Block[{Fgas,Fsolid},
Fgas[T_]:=3-1.5T;
Fsolid[T_]:=2-0.5T;
Graphics[{AbsoluteThickness[1],{Arrowheads[0.05],{Arrow[{{0,0},#1}],Text["\!\( "<>#2<>"\)",#1,#3]}&@@@{{{2.2,0},"T",{1,1.2}},{{0,3.2},"F",{2,1}}}},
First@Plot[{Fsolid[T],Fgas[T]},{T,0.2,1},PlotStyle->{Blue,Directive[Red,Dashed]}],
First@Plot[{Fsolid[T],Fgas[T]},{T,1,1.8},PlotStyle->{Directive[Blue,Dashed],Red}],
{Blue,Text["\!\(F\_solid\)",{0.2,Fsolid[0.2]},{-0.5,1.2}]},{Red,Text["\!\(F\_gas\)",{1.8,Fgas[1.8]},{-1.2,-0.2}]},
{Dotted,Line[{1,#}&/@{0,1.5}],Text["\!\(T\_*\)",{1,0},{0,1.2}]}},AspectRatio->0.5,ImageSize->180]]

The transition is generally a first order transition, because the first order derivative of the free energy ∂F/∂T is discontinuous at the transition.

Note: this simple argument assumes that each atom fluctuates in a mean-field potential, which is not valid in 1D and 2D. In low dimensions, collective fluctuations of atoms will destroy the crystal order (to be discussed later).

Lattice Structure

Lattice

A simple math model of lattice is a set of points as integer combinations of linearly independent primitive lattice vectors (Subscript[a, i]).

1D (n∈ℤ)

(Subscript[R, n])=n a.

2D ((Subscript[n, 1]),(Subscript[n, 2])∈ℤ)

(Subscript[R, n])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2]).

3D ((Subscript[n, 1]),(Subscript[n, 2]),(Subscript[n, 3])∈ℤ)

(Subscript[R, n])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])+(Subscript[n, 3])(Subscript[a, 3]).

Bravais lattice in dimension D is described by

Closure: if (Subscript[R, n]) and (Subscript[R, m]) are two sites on the Bravais lattice, (Subscript[R, n])+(Subscript[R, m])=(Subscript[R, n+m]) is also a site.

Homogeneity: every site on the Bravais lattice is equivalent. We can start with a representative site (say (Subscript[R, 0])=0) and reach all the other sites by translation.

General lattice (inhomogeneous): if the repeating structure contains more than one site, different types of sites (labeled by s) are further translated by different sub-lattice vectors (Subscript[r, s]).

Chain lattice (1D) (n∈ℤ)

a=1.

(Subscript[R, n])=n a=n.

Graphics[{Table[Translate[{{LightGray,AbsoluteThickness[1],Line[{{-0.5,0},{0.5,0}}]},{Black,AbsolutePointSize[5],Point[{0,0}]}},{n,0}],{n,-2,2}],
{Arrowheads[0.06],Arrow[{{0,0},{1,0}}]},Text["\!\(,a\)",{0.5,0},{0,-1}]},Axes->{True,False},AxesOrigin->{0,-1},Ticks->{Range[-2,2],None},
AspectRatio->0.2,ImageSize->180]

Square lattice (2D) ((Subscript[n, 1]),(Subscript[n, 2])∈ℤ)

({{{, ({{(Subscript[a, 1])=(1,0)}, {(Subscript[a, 2])=(0,1)}})}}).

(Subscript[R, n])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])=((Subscript[n, 1]),(Subscript[n, 2])).

Block[{L=2,a1,a2,Rs},
a1={1,0};
a2={0,1};
Rs=Tuples[Range[-L-1,L+1],2].{a1,a2};
Graphics[{{LightGray,AbsoluteThickness[1],Translate[Line[{{0,0},#}]&/@{a1,a2},#]&/@Rs},{Black,AbsolutePointSize[5],Point[Rs]},{Arrowheads[0.06],{Arrow[{{0,0},#1}],Text[StringTemplate["\!\(\*StyleBox[\"a\",Bold]\_`1`\)"][#2],#1,#3]}&@@@{{a1,"1",{0,0.8}},{a2,"2",{-1.4,0.2}}}}},PlotRange->L,PlotRangeClipping->True,PlotRangePadding->0.5,Frame->True,FrameTicks->{{Range[-L,L],{#,""}&/@Range[-L,L]},{Range[-L,L],{#,""}&/@Range[-L,L]}},ImagePadding->{{Automatic, 5}, {Automatic, 5}},ImageSize->160]]

Triangular lattice (2D) ((Subscript[n, 1]),(Subscript[n, 2])∈ℤ)

({{{, ({{(Subscript[a, 1])=(1,0)}, {(Subscript[a, 2])=((1/2),((Sqrt[3])/2))}})}}).

(Subscript[R, n])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])=((Subscript[n, 1])+((Subscript[n, 2])/2),(((Subscript[n, 2])(Sqrt[3]))/2)).

Block[{L=2,a1,a2,Rs},
a1={1,0};
a2={1,Sqrt[3]}/2;
Rs=Tuples[Range[-2L,2L],2].{a1,a2};
Graphics[{{LightGray,AbsoluteThickness[1],Translate[Line[{{0,0},#}]&/@CirclePoints[{1,0},3],#]&/@Rs},{Black,AbsolutePointSize[5],Point[Rs]},{Arrowheads[0.06],{Arrow[{{0,0},#1}],Text[StringTemplate["\!\(\*StyleBox[\"a\",Bold]\_`1`\)"][#2],#1,#3]}&@@@{{a1,"1",{0,0.8}},{a2,"2",{-1.4,0.2}}}}},PlotRange->L,PlotRangeClipping->True,PlotRangePadding->0.5,Frame->True,FrameTicks->{{Range[-L,L],{#,""}&/@Range[-L,L]},{Range[-L,L],{#,""}&/@Range[-L,L]}},ImagePadding->{{Automatic, 5}, {Automatic, 5}},ImageSize->160]]

Honeycomb lattice (2D) ((Subscript[n, 1]),(Subscript[n, 2])∈ℤ)

({{{, ({{(Subscript[a, 1])=((Sqrt[3]),0)}, {(Subscript[a, 2])=(((Sqrt[3])/2),(3/2))}})}}),({{{, ({{(Subscript[r, A])=((Sqrt[3]),1)}, {(Subscript[r, B])=(((Sqrt[3])/2),(1/2))}})}}).

(Subscript[R, n,A])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])+(Subscript[r, A])
=((Sqrt[3])(Subscript[n, 1])+(((Sqrt[3])(Subscript[n, 2]))/2)+(Sqrt[3]),((3(Subscript[n, 2]))/2)+1),
(Subscript[R, n,B])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])+(Subscript[r, B])
=((Sqrt[3])(Subscript[n, 1])+(((Sqrt[3])(Subscript[n, 2]))/2)+((Sqrt[3])/2),((3(Subscript[n, 2]))/2)+(1/2)).

Block[{L=1,a1,a2,Rs},
a1={Sqrt[3],0};
a2={Sqrt[3],3}/2;
Rs=Tuples[Range[-L-2,L+2],2].{a1,a2};
Graphics[{{LightGray,AbsoluteThickness[1],Translate[Line[{{0,0},#}]&/@CirclePoints[{1,π/2},3],#+{Sqrt[3],1}]&/@Rs},{AbsolutePointSize[5],Blue,Point[#+{Sqrt[3]/2,1/2}&/@Rs],Red,Point[#+{Sqrt[3],1}&/@Rs]},{Arrowheads[0.06],{Arrow[{{0,0},#1}],Text[StringTemplate["\!\(\*StyleBox[\"a\",Bold]\_`1`\)"][#2],#1,#3]}&@@@{{a1,"1",{0,0.8}},{a2,"2",{-1.,-0.2}}}},
{Blue,Text["\!\(,B\)",{Sqrt[3]/2,1/2},{1.5,0.5}],Red,Text["\!\(,A\)",{Sqrt[3],1},{-1.5,-0.5}]}},PlotRange->2L,PlotRangeClipping->True,PlotRangePadding->0.8,Frame->True,FrameTicks->{{Range[-2L,2L],{#,""}&/@Range[-2L,2L]},{Range[-2L,2L],{#,""}&/@Range[-2L,2L]}},ImagePadding->{{Automatic, 5}, {Automatic, 5}},ImageSize->160]]

Cubic lattice (3D) ((Subscript[n, 1]),(Subscript[n, 2]),(Subscript[n, 3])∈ℤ)

({{{, ({{(Subscript[a, 1])=(1,0,0)}, {(Subscript[a, 2])=(0,1,0)}, {(Subscript[a, 3])=(0,0,1)}})}}).

(Subscript[R, n])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])+(Subscript[n, 3])(Subscript[a, 3])=((Subscript[n, 1]),(Subscript[n, 2]),(Subscript[n, 3])).

Block[{L=1,a1,a2,a3,Rs},
a1={1,0,0};
a2={0,1,0};
a3={0,0,1};
Rs=Tuples[Range[-2L,2L],3].{a1,a2,a3};
Graphics3D[{Translate[{{Black,Sphere[{0,0,0},0.08]},{Black,Opacity[0.2],Tube[{{0,0,0},#},0.01]&/@{a1,a2,a3}}},#]&/@Rs,
{Black,Arrowheads[0.05],Arrow[Tube[{{0,0,0},#},0.02]]&/@{a1,a2,a3}},{Text[StringTemplate["\!\(\*StyleBox[\"a\",Bold]\_`1`\)"][#2],#1,#3]&@@@{{a1,1,{-1,1}},{a2,2,{1,-1}},{a3,3,{-1.2,1}}}}},PlotRange->L,PlotRangePadding->0.5,Axes->True,Ticks->{Range[-L,L],Range[-L,L],Range[-L,L]},Lighting->"Neutral",ViewPoint->{4,1,1},ImageSize->180]]

Diamond lattice (3D) ((Subscript[n, 1]),(Subscript[n, 2]),(Subscript[n, 3])∈ℤ)

({{{, ({{(Subscript[a, 1])=(1,1,0)}, {(Subscript[a, 2])=(0,1,1)}, {(Subscript[a, 3])=(1,0,1)}})}}),({{{, ({{(Subscript[r, A])=(0,0,0)}, {(Subscript[r, B])=((1/2),(1/2),(1/2))}})}}).

(Subscript[R, n,A])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])+(Subscript[n, 3])(Subscript[a, 3])+(Subscript[r, A])
=((Subscript[n, 1])+(Subscript[n, 3]),(Subscript[n, 1])+(Subscript[n, 2]),(Subscript[n, 2])+(Subscript[n, 3])),
(Subscript[R, n,B])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])+(Subscript[n, 3])(Subscript[a, 3])+(Subscript[r, B])
=((Subscript[n, 1])+(Subscript[n, 3])+(1/2),(Subscript[n, 1])+(Subscript[n, 2])+(1/2),(Subscript[n, 2])+(Subscript[n, 3])+(1/2)).

Block[{L=1,a1,a2,a3,RAs,RBs},
a1={1,1,0};
a2={0,1,1};
a3={1,0,1};
RAs=Tuples[Range[-2L,2L],3].{a1,a2,a3};
RBs=#+{1,1,1}/2&/@RAs;
Graphics3D[{{Red,Sphere[#,0.1]&/@RAs,Blue,Sphere[#,0.1]&/@RBs},{Gray,Translate[Tube[{{0,0,0},#}&/@{{1,-1,1},{-1,1,1},{1,1,-1},{-1,-1,-1}}/2,0.015],#]&/@RBs},
{Black,Arrowheads[0.06],Arrow[Tube[{{0,0,0},#},0.02]]&/@{a1,a2,a3}},
{Text[StringTemplate["\!\(\*StyleBox[\"a\",Bold]\_`1`\)"][#2],#1,#3]&@@@{{a1,1,{0,1}},{a2,2,{1.5,-0.2}},{a3,3,{-1,1}}}},{Translate[GraphicsComplex[{{0,0,0},a1,a2,a3},{Opacity[0.2],EdgeForm[None],Polygon[{{1,2,3},{1,3,4},{1,2,4},{2,3,4}}]}],#]&/@{{0,0,0},-a1,-a2,-a3}},
{Red,Text["\!\(,A\)",{0,0,0},{0.2,-1}],Blue,Text["\!\(,B\)",{1,1,1}/2,{0,1.2}]}},PlotRange->L,PlotRangePadding->0.2,Axes->True,Ticks->{Range[-L,L],Range[-L,L],Range[-L,L]},Lighting->"Neutral",ViewPoint->{2 ,-5,1},ImageSize->180]]

Unit Cell

Block[{Rs},
Rs=Tuples[{Range[-5,5],Range[-3,3]}].{{1,0},{1,Sqrt[3]}/2};
Graphics[{{LightGray,Parallelogram[{-5/2,-Sqrt[3]/2},{{1,0},{0,Sqrt[3]}}],Parallelogram[{-1,0},{{1,0},{1/2,Sqrt[3]/2}}],Parallelogram[{3/2,Sqrt[3]/2},{{1,0},{-1/2,Sqrt[3]/2}}],Translate[Polygon[CirclePoints[{1/Sqrt[3],π/6},6]],{3/2,-Sqrt[3]/2}]},{Black,AbsolutePointSize[5],Point[Rs]}},PlotRange->{{-3.75,3.75},{-2,2}}]]

Primitive unit cell: a unit cell of minimal volume, e.g. parallelogram (parallelepiped) spanned by primitive lattice vectors (Subscript[a, i]).

Block[{Rs,puc},
Rs=Tuples[{Range[-5,5],Range[-3,3]}].{{1,0},{1,Sqrt[3]}/2};
puc[p_,vs_]:={Parallelogram[p,vs],{Black,Translate[Arrow[{{0,0},#}&/@vs],p]}};
Graphics[{{LightGray,puc[{-1,0},{{1,0},{1/2,Sqrt[3]/2}}],puc[{3/2,Sqrt[3]/2},{{1,0},{-1/2,Sqrt[3]/2}}],Translate[Polygon[CirclePoints[{1/Sqrt[3],π/6},6]],{3/2,-Sqrt[3]/2}]},{Black,AbsolutePointSize[5],Point[Rs]}},PlotRange->{{-3.75,3.75},{-2,2}}]]

The choice of primitive lattice vectors for a lattice is not unique, hence the primitive unit cell is also not unique.

However, the volume of the primitive unit cell is unique, which is given by the absolute value of the determinant of the matrix A formed by all primitive lattice vectors (see Wikipedia Parallelepiped:Volume)

In 2D, Ω=|(Subscript[a, 1])×(Subscript[a, 2])|.

In 3D, Ω=|(Subscript[a, 1])·((Subscript[a, 2])×(Subscript[a, 3]))|.

Calculate the volume of primitive unit cell for the triangular lattice defined in Eq. ().

({{{, ({{(Subscript[a, 1])=(1,0)}, {(Subscript[a, 2])=((1/2),((Sqrt[3])/2))}})}})=>(({{(Subscript[a, 1])}, {(Subscript[a, 2])}}))=(({{1, 0}, {(1/2), ((Sqrt[3])/2)}})),

Ω=det(({{1, 0}, {(1/2), ((Sqrt[3])/2)}}))=((Sqrt[3])/2).

Wigner-Seitz cell (Voronoi cell): a primitive unit cell enclosed by perpendicular bisectors between a Bravais lattice point and all its neighbors.

Block[{Rs,puc},
Rs=Tuples[{Range[-5,5],Range[-3,3]}].{{1,0},{1,Sqrt[3]}/2};
Graphics[{{Translate[{LightGray,Polygon[CirclePoints[{1/Sqrt[3],π/6},6]],{Black,Dashing[{0.015,0.015},0.002],Line[{{0,0},#}&/@CirclePoints[6]]},{Black,AbsoluteThickness[1],Rotate[Line[{{1/2,-1/2},{1/2,1/2}}],# π/3,{0,0}]&/@Range[6]}},{3/2,-Sqrt[3]/2}]},,{Black,AbsolutePointSize[5],Point[Rs]}},PlotRange->{{-3.75,3.75},{-2,2}}]]

Lattice Symmetry

General form of coordinate transformation (affine transformation)

R→R'=R M+b,

where

M - matrix (to implement rotations, reflections, inversions),

b - vector (to implement translation).

Translation Symmetry (M=𝟙)

(Subscript[R, n,s])→(Subscript[R, n,s]^′)=(Subscript[R, n,s])+(Subscript[R, m])=(Subscript[R, n+m,s])

Point Group Symmetry (b=0)

Rotation by angle θ along an axis. The rotation is said to be k-fold if the angle is θ=2π/k. In 2D, rotation axis reduces to a point, and the rotation about the origin is implemented by

R→R(({{cos θ, sin θ}, {-sin θ, cos θ}})).

In 3D, a rotation about axis N (specified by a unit vector) by θ angle is implemented by

R→R cos θ+(R·N)N(1-cos θ)+N×R sin θ.

Reflection about a mirror plane (in 3D; reduces to a line in 2D). The mirror plane can be specified by the normal vector N (the unit vector perpendicular to the mirror plane).

R→R-2(R·N)N.

Show the composition of two reflections induced by normal vectors (Subscript[N, 1]) and (Subscript[N, 2]) is a rotation along the axis ((Subscript[N, 1])×(Subscript[N, 2]))/sin θ by the angle 2θ, where θ is the angle between (Subscript[N, 1]) and (Subscript[N, 2]) (i.e. (Subscript[N, 1])·(Subscript[N, 2])=cos θ).

R→R-2(R·(Subscript[N, 1]))(Subscript[N, 1])
→(R-2(R·(Subscript[N, 1]))(Subscript[N, 1]))-2((R-2(R·(Subscript[N, 1]))(Subscript[N, 1]))·(Subscript[N, 2]))(Subscript[N, 2])
=R-2(R·(Subscript[N, 1]))(Subscript[N, 1])-2(R·(Subscript[N, 2]))(Subscript[N, 2])+4(R·(Subscript[N, 1]))((Subscript[N, 1])·(Subscript[N, 2]))(Subscript[N, 2]).

R→R cos 2θ+R·((Subscript[N, 1])×(Subscript[N, 2]))((Subscript[N, 1])×(Subscript[N, 2]))((1-cos 2θ)/((sin^2) θ))+((Subscript[N, 1])×(Subscript[N, 2]))×R ((sin 2θ)/(sin θ))
=R (2 (cos^2)θ-1)+2R·((Subscript[N, 1])×(Subscript[N, 2]))((Subscript[N, 1])×(Subscript[N, 2]))+2R×((Subscript[N, 1])×(Subscript[N, 2])) cos θ
=R (2 ((((Subscript[N, 1])·(Subscript[N, 2])))^2)-1)+2R·((Subscript[N, 1])×(Subscript[N, 2]))((Subscript[N, 1])×(Subscript[N, 2]))+2 ((Subscript[N, 1])·(Subscript[N, 2]))((Subscript[N, 1])×(Subscript[N, 2]))×R

R=(Subscript[R, ||])+(Subscript[R, ⊥]).

(Subscript[R, ⊥])=(1/((sin^2)θ))R·((Subscript[N, 1])×(Subscript[N, 2]))((Subscript[N, 1])×(Subscript[N, 2]))
=((R·((Subscript[N, 1])×(Subscript[N, 2]))((Subscript[N, 1])×(Subscript[N, 2])))/(1-((((Subscript[N, 1])·(Subscript[N, 2])))^2))).

(Subscript[R, ||])=(Subscript[x, 1])(Subscript[N, 1])+(Subscript[x, 2])(Subscript[N, 2]),

R·(Subscript[N, 1])=(Subscript[R, ||])·(Subscript[N, 1])=(Subscript[x, 1])+(Subscript[x, 2])((Subscript[N, 1])·(Subscript[N, 2])),
R·(Subscript[N, 2])=(Subscript[R, ||])·(Subscript[N, 2])=(Subscript[x, 2])+(Subscript[x, 1])((Subscript[N, 1])·(Subscript[N, 2])),

(Subscript[x, 1])=((R·(Subscript[N, 1])-((Subscript[N, 1])·(Subscript[N, 2]))(R·(Subscript[N, 2])))/(1-((((Subscript[N, 1])·(Subscript[N, 2])))^2))),
(Subscript[x, 2])=((R·(Subscript[N, 2])-((Subscript[N, 1])·(Subscript[N, 2]))(R·(Subscript[N, 1])))/(1-((((Subscript[N, 1])·(Subscript[N, 2])))^2))),

(Subscript[R, ||])=(1/(1-((((Subscript[N, 1])·(Subscript[N, 2])))^2)))((R·(Subscript[N, 1])-((Subscript[N, 1])·(Subscript[N, 2]))(R·(Subscript[N, 2])))(Subscript[N, 1])+(R·(Subscript[N, 2])-((Subscript[N, 1])·(Subscript[N, 2]))(R·(Subscript[N, 1])))(Subscript[N, 2])).

R(1-((((Subscript[N, 1])·(Subscript[N, 2])))^2))=(R·(Subscript[N, 1])-((Subscript[N, 1])·(Subscript[N, 2]))(R·(Subscript[N, 2])))(Subscript[N, 1])+(R·(Subscript[N, 2])-((Subscript[N, 1])·(Subscript[N, 2]))(R·(Subscript[N, 1])))(Subscript[N, 2])+R·((Subscript[N, 1])×(Subscript[N, 2]))((Subscript[N, 1])×(Subscript[N, 2])).

((Subscript[N, 1])×(Subscript[N, 2]))×R=(R·(Subscript[N, 1]))(Subscript[N, 2])-(R·(Subscript[N, 2]))(Subscript[N, 1]).

R→R (2 ((((Subscript[N, 1])·(Subscript[N, 2])))^2)-1)+2(R(1-((((Subscript[N, 1])·(Subscript[N, 2])))^2))-(R·(Subscript[N, 1])-((Subscript[N, 1])·(Subscript[N, 2]))(R·(Subscript[N, 2])))(Subscript[N, 1])-(R·(Subscript[N, 2])-((Subscript[N, 1])·(Subscript[N, 2]))(R·(Subscript[N, 1])))(Subscript[N, 2]))+2 ((Subscript[N, 1])·(Subscript[N, 2]))((R·(Subscript[N, 1]))(Subscript[N, 2])-(R·(Subscript[N, 2]))(Subscript[N, 1]))
= R-2(R·(Subscript[N, 1]))(Subscript[N, 1])-2(R·(Subscript[N, 2]))(Subscript[N, 2])+4 (R·(Subscript[N, 1]))((Subscript[N, 1])·(Subscript[N, 2]))(Subscript[N, 2]).

Inversion about an inversion center. In 2D an inversion is equivalent to a 2-fold rotation. Inversion about the origin is implemented by

R→-R.

If the inversion center (Subscript[r, c]) is not at the origin, how does the lattice coordinate transform under the inversion?

R(Overscript[→, translate])R-(Subscript[r, c])(Overscript[→, inverse])(Subscript[r, c])-R(Overscript[→, (translate^(-1))])2(Subscript[r, c])-R.

Nonsymmorphic Symmetry

Glide reflection: translate along an axis followed by a reflection about the same axis.

Example: Shastry-Sutherland lattice.

Manipulate[Block[{sub,pts,bds},
sub=AssociationThread[Range@Length@#->#&[{{1-δ,1-δ},{1+δ,3-δ},{3-δ,1+δ},{3+δ,3+δ}}/4]];
pts=Plus@@@Tuples[{Tuples[Range[-1,1],2],Values@sub}];
bds=
GroupBy[ReplaceList[pts,{___,p1:{_,_},___,p2:{_,_},___}/;Norm[p1-p2]≤1/Sqrt[2]:>{p1,p2}],Abs[Norm[#[[1]]-#[[2]]]-(1-δ)/Sqrt[2]]<0.001&];
Graphics[{
{Table[{AbsoluteThickness@If[typ,2,1.5],If[typ,Gray,LightGray],Line[#]&/@bds[typ]},{typ,{True,False}}]},{Black,AbsolutePointSize[5],Point[pts]},{AbsoluteThickness[1],Arrowheads[0.05],Arrow[{#,1/4}&/@{-1,2}]}},ImageSize->{Automatic,140}]
],{{δ,0.5},0.001,0.5,ImageSize->Small},Paneled->False,SaveDefinitions->True]

Screw displacement: translation along an axis followed by a rotation around the same axis.

Block[{pts},
pts=Table[{x,Cos[π x/6],Sin[π x/6]},{x,-6,6}];
Graphics3D[{Rotate[{#2,Sphere[#,0.2]&/@pts,Tube[pts,0.1]},#1,{1,0,0},{0,0,0}]&@@@{{0,Red},{π,Blue}},
Line[{#,{1,-1,-1}#}&/@pts],{Arrow[{{-7,0,0},{7,0,0}}]}},ViewPoint->{0,-2,0.5},ImageSize->240]]

Reciprocal Lattice

Wave Propagating on Lattice

(Subscript[ψ, k])(r)=(e^(i k·r))

length of k - wave number |k|=2π/λ, reciprocal to the wave length λ,

direction of k - wave propagation direction.

Block[{L=2.2,Rs,b1,b2,plt,ms},
Rs=Cases[Tuples[Range[-3,3],2].{{1,0},{1/2,Sqrt[3]/2}},{x_,y_}/;Abs[x]<=2&&Abs[y]<=2];
b1={0,4π/Sqrt[3]};
b2={-2π,2π/Sqrt[3]};
plt[fun_]:=Graphics[{Raster[Table[fun[{x,y}],{y,-L,L,0.1},{x,-L,L,0.1}],{{-L,-L},{L,L}},ColorFunction->ComplexColor]},BaselinePosition->(Axis->Baseline),ImageSize->80];
ms={{0,1/3},{0,4/3},{1,1/3},{1,-2/3}};
Row[Table[
Block[{k=m.{b1,b2}},Show[{plt[Exp[I k.#]&],Graphics[{Arrowheads[0.12],If[Norm@k>0,Arrow[{0.1#,1#}&@Normalize[k]],Nothing]}]},PlotLabel->Style[MatrixForm[{k}],12,Black,FontFamily->"LaTeX"]]],{m,ms}]]
]

Block[{L=2.2,Rs,b1,b2,plt,ms},
Rs=Cases[Tuples[Range[-3,3],2].{{1,0},{1/2,Sqrt[3]/2}},{x_,y_}/;Abs[x]<=2&&Abs[y]<=2];
b1={0,4π/Sqrt[3]};
b2={-2π,2π/Sqrt[3]};
plt[fun_]:=Graphics[{Raster[Table[fun[{x,y}],{y,-L,L,0.1},{x,-L,L,0.1}],{{-L,-L},{L,L}},ColorFunction->ComplexColor],
{AbsoluteThickness[1],Circle[#,Offset@{2,2}]&/@Rs}},BaselinePosition->(Axis->Baseline),ImageSize->80];
ms={{0,1/3},{0,4/3},{1,1/3},{1,-2/3}};
Row[Table[
Block[{k=m.{b1,b2}},Show[{plt[Exp[I k.#]&],Graphics[{Arrowheads[0.12],If[Norm@k>0,Arrow[{0.1#,1#}&@Normalize[k]],Nothing]}]},PlotLabel->Style[MatrixForm[{k}],12,Black,FontFamily->"LaTeX"]]],{m,ms}]]
]

Block[{L=2.2,Rs,b1,b2},
Rs=Cases[Tuples[Range[-3,3],2].{{1,0},{1/2,Sqrt[3]/2}},{x_,y_}/;Abs[x]<=2&&Abs[y]<=2];
b1={0,4π/Sqrt[3]};
b2={-2π,2π/Sqrt[3]};
Block[{k={0,1/3}.{b1,b2}},Graphics[{
{EdgeForm[AbsoluteThickness[1]],{ComplexColor[Exp[I k.#]],Disk[#,Offset@{2,2}]}&/@Rs}},BaselinePosition->(Axis->Baseline),ImageSize->80]]
]

Two waves of wave vectors k and k' are equivalent (indistinguishable) on a lattice, if their wave momentum difference G=k'-k satisfies

(e^(i G·R))=1

for all lattice points R.

Because, in this case,

(Subscript[ψ, k'])(R)=(e^(i k'·R))=(e^(i (k+G)·R))=(e^(i k·R))(e^(i G·R))=(Subscript[ψ, k])(R).

Note: this does not mean (Subscript[ψ, k'])(r)=(Subscript[ψ, k])(r) for any point r in the space, but the equality only holds for discrete lattice points R.

The solutions G of Eq. () are called lattice momenta of the lattice. It turns out that all solutions of G also forms a lattice in the momentum space, called the reciprocal lattice.

Reciprocal Lattice

Recall Eq. (): every Bravais lattice vector (Subscript[R, n])=(∑_(i=1)^D)(Subscript[n, i])(Subscript[a, i]). is an integer combination of primitive lattice vectors (Subscript[a, i]). Then (e^(i G·(Subscript[R, n])))=1 implies

G·(Subscript[R, n])=(∑_(i=1)^D)(Subscript[n, i])G·(Subscript[a, i])∈2π ℤ,

and this must hold for any choice of n.

Consider special cases, when n is a unit vector

n=(1,0,…)=>G·(Subscript[a, 1])∈2π ℤ,
n=(0,1,…)=>G·(Subscript[a, 2])∈2π ℤ,
… (enumerating over all dimensions).

In summary,

G·(Subscript[a, j])=2π (Subscript[m, j])

with (Subscript[m, j])∈ℤ (for j=1,2,…,D). This is a necessary and sufficient condition for G to be a reciprocal lattice point

Necessity: if any condition in Eq. () is not satisfied, we can take the corresponding choice of n to construct a violation of (e^(i G·(Subscript[R, n])))=1.

Sufficiency: as long as all conditions in Eq. () is satisfied, we have

G·(Subscript[R, n])=2π ((∑_(i=1)^D)(Subscript[m, i])(Subscript[n, i]))∈2π ℤ,

such that (e^(i G·(Subscript[R, n])))=1.

Now we define a matrix A by collecting rows of primitive lattice vectors (Subscript[a, i])

A(Overscript[=, def])(({{(Subscript[a, 1])}, {(Subscript[a, 2])}, {…}})),

Eq. () can then be written in the linear algebra form

G Aᵀ=2π m.

G=((Subscript[G, 1]),(Subscript[G, 2]),…)∈(ℝ^D) - a solution of reciprocal vector

Aᵀ - transpose of A matrix, such that each row vector of A (other than column vector) is dot producing with the G vector.

m=((Subscript[m, 1]),(Subscript[m, 2]),…)∈(ℤ^D) - uniquely labels one solution of (Subscript[G, m]) via Eq. (). Components of m are also called Miller indices.

The solution of (Subscript[G, m]) can be formally obtained by matrix inversion.

(Subscript[G, m])=2π m (((Aᵀ))^(-1))=m B,

where we introduce the matrix B of

B(Overscript[=, def])(({{(Subscript[b, 1])}, {(Subscript[b, 2])}, {…}}))=2π (((Aᵀ))^(-1)).

Eq. () implies that the solution (Subscript[G, m]) is given by

which also takes the form of a lattice => the reciprocal lattice.

These (Subscript[b, j]) vectors actually form the primitive lattice vectors of the reciprocal lattice.

Eq. () indicates that the primitive vectors of the real and reciprocal lattices are related by B Aᵀ=2π 𝟙, or more explicitly as

The name “reciprocal” becomes clear now, as their primitive lattice vectors are reciprocal (as matrix inversion) to each other. When the direct lattice expands, the reciprocal lattice will shrink, and vice versa.

1D lattice

2D lattice Ω=|(Subscript[a, 1])×(Subscript[a, 2])|,

where (OverHat[z]) is the unit vector normal to the 2D lattice plane.

3D lattice Ω=|(Subscript[a, 1])·((Subscript[a, 2])×(Subscript[a, 3]))|,

Find the primitive lattice vectors of the reciprocal lattice of the triangular lattice as defined by Eq. (). Plot the reciprocal lattice.

({{{, ({{(Subscript[a, 1])=(1,0)}, {(Subscript[a, 2])=((1/2),((Sqrt[3])/2))}})}}),

(Subscript[b, 1])=((2π)/Ω) (Subscript[a, 2])×(OverHat[z])
=((4π)/(Sqrt[3]))((1/2),((Sqrt[3])/2),0)×(0,0,1)
=((4π)/(Sqrt[3]))(((Sqrt[3])/2),-(1/2),0)
=(2π,-((2π)/(Sqrt[3]))).

(Subscript[b, 2])=-((2π)/Ω) (Subscript[a, 1])×(OverHat[z])
=-((4π)/(Sqrt[3]))(1,0,0)×(0,0,1)
=((4π)/(Sqrt[3]))(0,1,0)
=(0,((4π)/(Sqrt[3]))),

Block[{b1,b2,Gs,Q=4π},
b1={2π,-2π/Sqrt[3]};
b2={0,4π/Sqrt[3]};
Gs=Tuples[Range[-3,3],2].{b1,b2};
Graphics[{Red,Point[Gs]},PlotRange->Q,PlotRangePadding->Scaled[.02],PlotRangeClipping->True,Frame->True,FrameStyle->Black,FrameTicks->{{Range[-Q,Q,2π],None},{Range[-Q,Q,2π],None}},ImageSize->{Automatic,120}]]

Block[{L=12,Q=4π,Rs,rs,qs},
Rs=Cases[Tuples[Range[-14,14],2].{{1,0},{1/2,Sqrt[3]/2}},p_/;Norm[p]<=L];
qs=N@Table[{qx,qy},{qy,Q,-Q,-2π/L},{qx,-Q,Q,2π/L}];
ComplexMatrixPlot[Total[Exp[I tTr[{qs,Rs},{{3,5}}]],{-1}]/Length[Rs],DataRange->{{-1,1},{-1,1}}Q,Frame->True,FrameStyle->Black,FrameTicks->{Range[-Q,Q,2π],Range[-Q,Q,2π]},PlotLabel->Style["\!\(ρ\&~(q)\)",14,Black],ImageSize->{Automatic,140}]]

Brillouin Zone

The first Brillouin zone (1st BZ): the Wigner-Seitz cell centered at G=0 on the reciprocal lattice; also as the set of all points k in the momentum space, that are closest to G=0 than another points G on the reciprocal lattice.

The nth Brillouin zone: the set of all points k in the momentum space, that are nth closest to G=0 among another points G on the reciprocal lattice.

On 1D reciprocal lattice

Block[{M=2,h=0.4},
Graphics[{{Table[Block[{col=ColorData[0][2m],BZ=<|1->"1st",2->"2nd",3->"3rd",4->"4th",5->"5th"|>[2m]},{{Lighter[col,0.7],Rectangle[{-m,-h},{m,h}]},{col,If[m==1/2,Text[BZ<>" BZ",{0,h},{0,1}],{Text[BZ,{m-1/4,h},{0,1}],Text[BZ,{-m+1/4,h},{0,1}]}]}}],{m,M+1/2,1/2,-1/2}]},
{Arrow[{#(M+0.5),0}&/@{-1,1}]},
{AbsolutePointSize[5],Table[Translate[{Point[{0,0}],Text[StringTemplate["\!\(G\_`1`\)"][m],{0,0},{0,1.2}]},{m,0}],{m,-M,M}]},
{AbsoluteThickness[1],Opacity[0.5],Table[If[m≠0,Translate[{Line[{0,# h}&/@{-1,1}]},{m,0}]],{m,-M,M,1/2}]},Text["\!\( k\)",{M+1/2,0},{-1.2,0}]},ImageSize->300]]

On 2D square reciprocal lattice

Block[{M=4,pts,k2p,polyreg,update},
pts=Tuples[Range[-M,M],2];
k2p=Block[{n2p=GroupBy[pts,Norm]},n2p/@AssociationThread[Range[0,Length[#]-1]->#&@SortBy[Keys@n2p,N]]];
polyreg[poly_]:=ReplaceRepeated[poly,{h:(Polygon|Triangle)[{pre___,p1_,p2_,post___}]/;Norm[p1-p2]<0.001:>Head[h][{pre,(p1+p2)/2,post}],h:(Polygon|Triangle)[{p1_,mid___,p2_}]/;Norm[p1-p2]<0.001:>Head[h][{(p1+p2)/2,mid}]}];
update[rank_,p_List]:=Block[{inner,outer,ins,mids,outs},
{inner,outer}=HalfPlane[(0.5-0.00001#)p,{{0,-1},{1,0}}.p,# p]&/@{-1,1};
{ins,mids,outs}=Flatten[{{Lookup[#,True,{}]},{Lookup[#,False,{}],Lookup[#,True,{}]}&@GroupBy[Lookup[#,False,{}],RegionWithin[outer,#]&]}&@GroupBy[Keys@rank,RegionWithin[inner,#]&],1];
Association[Join[Table[poly->rank[poly],{poly,ins}],Table[polyreg@RegionIntersection[poly,inner]->rank[poly],{poly,mids}],Table[polyreg@RegionIntersection[poly,outer]->rank[poly]+1,{poly,mids}],Table[poly->rank[poly]+1,{poly,outs}]]]];
update[rank_,k_Integer]:=Fold[update,rank,k2p@k];
Block[{rank},
rank=Fold[update,<|Polygon[CirclePoints[{1.7,0},8]]->1|>,Range[6]];
rank=Select[rank,#<=8&];
Graphics[{{EdgeForm[Directive[AbsoluteThickness[1],Opacity[0.2]]],Table[Block[{col={Blue,Red,Green,Yellow,Purple,Cyan,Orange,LightGreen}[[rank[poly]]]},{{Lighter[col,0.6],poly}(*,If[Area[poly]>=1/8,{col,Text[closerank[poly],RegionCentroid[poly]]}]*)}],{poly,Keys@rank}]},
{Point[pts]}},PlotRange->2,PlotRangePadding->0.2,ImageSize->160]
]]

On 2D triangular reciprocal lattice

Block[{M=3,pts,k2p,polyreg,update},
pts=N@Cases[Tuples[Range[-M,M],2].{{Sqrt[3]/2,-1/2},{0,1}},p_/;Norm[p]<M];
k2p=Block[{n2p=GroupBy[pts,Norm]},n2p/@AssociationThread[Range[0,Length[#]-1]->#&@SortBy[Keys@n2p,N]]];
polyreg[poly_]:=ReplaceRepeated[poly,{h:(Polygon|Triangle)[{pre___,p1_,p2_,post___}]/;Norm[p1-p2]<0.001:>Head[h][{pre,(p1+p2)/2,post}],h:(Polygon|Triangle)[{p1_,mid___,p2_}]/;Norm[p1-p2]<0.001:>Head[h][{(p1+p2)/2,mid}]}];
update[rank_,p_List]:=Block[{inner,outer,ins,mids,outs},
{inner,outer}=HalfPlane[(0.5-0.00001#)p,{{0,-1},{1,0}}.p,# p]&/@{-1,1};
{ins,mids,outs}=Flatten[{{Lookup[#,True,{}]},{Lookup[#,False,{}],Lookup[#,True,{}]}&@GroupBy[Lookup[#,False,{}],RegionWithin[outer,#]&]}&@GroupBy[Keys@rank,RegionWithin[inner,#]&],1];
Association[Join[Table[poly->rank[poly],{poly,ins}],Table[polyreg@RegionIntersection[poly,inner]->rank[poly],{poly,mids}],Table[polyreg@RegionIntersection[poly,outer]->rank[poly]+1,{poly,mids}],Table[poly->rank[poly]+1,{poly,outs}]]]];
update[rank_,k_Integer]:=Fold[update,rank,k2p@k];
Block[{rank},
rank=Fold[update,<|Polygon[CirclePoints[{Sqrt[3],0},6]]->1|>,Range[4]];
rank=Select[rank,#<=8&];
Graphics[{{EdgeForm[Directive[AbsoluteThickness[1],Opacity[0.2]]],Table[Block[{col={Blue,Red,Green,Yellow,Purple,Cyan,Orange,LightGreen}[[rank[poly]]]},{{Lighter[col,0.6],poly}(*,If[Area[poly]>=1/8,{col,Text[closerank[poly],RegionCentroid[poly]]}]*)}],{poly,Keys@rank}]},
{Point[pts]}},PlotRange->2,PlotRangePadding->0.2,ImageSize->160]
]]

The 1st BZ provides a compete set of inequivalent wave vectors. Any wave vector outside the 1st BZ will be equivalent to one inside the 1st BZ (related by certain lattice momentum).

Every BZ has exactly the same volume (in the momentum space):

det|B|=((((2π))^n)/(det|A|))=((((2π))^n)/Ω).

Higher BZ can be folded to the 1st BZ.

Any physical property (e.g. frequency, group velocity) of the wave propagating on a lattice can be fully described within the 1st BZ.

Smith, Philip. Geometric and Topological Methods for Applications to Materials and Data Skeletonisation. PhD thesis, University of Liverpool (2021).

Scattering Experiment

Elastic Scattering

X-ray scattering.

Neutron scattering.

Electron scattering.

Real Space Picture

Plane wave hits an atom, the atom responses with a spherical wave.

Block[{R=3,r=0.1,in,out,pts,plt},
in[{x_,y_}]:=Exp[2π I x];
out[{x_,y_}]:=HankelH1[0,2π Sqrt[x^2+y^2]];
pts={{0,0}};plt[fun_]:=Graphics[{Raster[Table[fun[{x,y}],{y,-R,R,0.1},{x,-R,R,0.1}],{{-R,-R},{R,R}},ColorFunction->ComplexColor],{EdgeForm[Directive[Black,AbsoluteThickness[1]]],{ComplexColor[fun[#+{r,0}]],Disk[#,r]}&/@pts}},ImageSize->90];
Style[Row[{plt[0.7in[#]&],"+",plt[-0.7Sum[in[pt]out[#-pt],{pt,pts}]&],"=",plt[0.7(in[#]-Sum[in[pt]out[#-pt],{pt,pts}])&]}],14,FontFamily->"LaTeX"]]

Each atom behaves like an independent scatterer for the incident wave (to the leading order, ignoring the fact that the atoms can continue to scatter the wave that was scattered from other atoms).

Block[{R=3,r=0.1,in,out,pts,plt},
in[{x_,y_}]:=Exp[2π I x];
out[{x_,y_}]:=HankelH1[0,2π Sqrt[x^2+y^2]];
pts={{0,1},{0,-1}};plt[fun_]:=Graphics[{Raster[Table[fun[{x,y}],{y,-R,R,0.1},{x,-R,R,0.1}],{{-R,-R},{R,R}},ColorFunction->ComplexColor],{EdgeForm[Directive[Black,AbsoluteThickness[1]]],{ComplexColor[fun[#+{r,0}]],Disk[#,r]}&/@pts}},ImageSize->90];
Style[Row[{plt[0.7in[#]&],"+",plt[-0.7Sum[in[pt]out[#-pt],{pt,pts}]&],"=",plt[0.7(in[#]-Sum[in[pt]out[#-pt],{pt,pts}])&]}],14,FontFamily->"LaTeX"]]

If there is an array of atoms, the scattering wave emitted from each atom can interfere → leading to Bragg diffraction.

Block[{R=3,r=0.1,in,out,pts,plt},
in[{x_,y_}]:=Exp[2π I x];
out[{x_,y_}]:=Exp[2π I Sqrt[x^2+y^2]]/((x^2+y^2)/r^2+1)^(1/4);
pts=Flatten[Table[{x,y},{y,-1,1},{x,-1,1}],1];plt[fun_]:=Graphics[{Raster[Table[fun[{x,y}],{y,-R,R,0.1},{x,-R,R,0.1}],{{-R,-R},{R,R}},ColorFunction->ComplexColor],{EdgeForm[Directive[Black,AbsoluteThickness[1]]],{ComplexColor[fun[#+{r,0}]],Disk[#,r]}&/@pts}},ImageSize->90];
Style[Row[{plt[0.5in[#]&],"+",plt[-0.5Sum[in[pt]out[#-pt],{pt,pts}]&],"=",plt[0.5(in[#]-Sum[in[pt]out[#-pt],{pt,pts}])&]}],14,FontFamily->"LaTeX"]]

As we zoom out, the wave scattered off the crystal looks like beams of plane waves along special directions.

Manipulate[Block[{R=8,r=0.1,in,out,pts,plt},
in[{x_,y_}]:=Exp[2π I x];
out[{x_,y_}]:=Exp[2π I Sqrt[x^2+y^2]]/((x^2+y^2)/r^2+1)^(1/4);
pts=Flatten[Table[ReIm[(x+I y)Exp[I θ]],{y,-1,1},{x,-1,1}],1];plt[fun_]:=Graphics[{Raster[Table[fun[{x,y}],{y,-R,R,0.2},{x,-R,R,0.2}],{{-R,-R},{R,R}},ColorFunction->ComplexColor],{EdgeForm[Directive[Black,AbsoluteThickness[1]]],{ComplexColor[fun[#]],Disk[#,r]}&/@pts}},ImageSize->150];
Show[plt[-0.7Sum[in[pt]out[#-pt],{pt,pts}]&],Graphics[{Arrowheads[0.06],Arrow[{2#,7#}]&/@ReIm[Cases[1+(Tuples[Range[-2,2],{2}].{1,I})Exp[I θ],z_/;Abs[Abs[z]-1]<0.15]]}]]],{θ,0,π/2,ImageSize->Small},ContinuousAction->False,Paneled->False,SaveDefinitions->True]

Momentum Space Picture

Plane wave of wave momentum k (length of k - 2π/λ, direction of k - propagation direction)

(Subscript[ψ, k])(r)=(e^(i k·r)).

Spherical wave = superposition of plane waves of the same wave length in all directions

(Subscript[ϕ, k])(r)=(∫_(|k|=k))(Subscript[ψ, k])(r)=(∫_(|k|=k))(e^(i k·r))
=(1/(2π))∫dθ (e^(i k r cos θ))=(Subscript[J, 0])(k r).

An atom scatters a plane wave of momentum k to a spherical wave of the same wave number k ⇔ to a equal-probability superposition of plane waves of all possible momentum k', as long as |k'|=|k|=k (the wave length does not change, as it is fixed by the wave frequency).

Graphics[{{Arrowheads[0.08],{Red,Arrow[{{-2,0},{-0.5,0}}],Text["\!\( k\)",{-2,0},{1.2,0}]},{Blue,Arrow[{0.5#,2#}]&/@CirclePoints[{1,π/12},12],Text["\!\( k\^′\)",{0,1.5}]}},{EdgeForm[AbsoluteThickness[1]],White,Disk[{0,0},0.1]}},ImageSize->100]

Γ(k→k')∝((|(Bra[k'])ρ(Ket[k])|)^2)δ((Subscript[E, k'])-(Subscript[E, k])),

the scattering amplitude (Bra[k'])ρ(Ket[k]) is given by the wave function overlap

(Bra[k'])ρ(Ket[k])=(∫_r)(Subscript[ψ, k']^( * ))(r)ρ(r)(Subscript[ψ, k])(r)
=(∫_r)ρ(r)(e^(-i(k'-k)·r))
=(OverTilde[ρ])(k'-k)

(Subscript[ψ, k])(r) - incident wave,

ρ(r) - scatterer density distribution (~ atomic density distribution),

ρ(r)(Subscript[ψ, k])(r) - incident wave enveloped by the atomic density (no atom no scattering, more atom stronger scattering),

(∫_r)(Subscript[ψ, k']^( * ))(r)… - inner product with (Subscript[ψ, k'])(r) to extract the amplitude of wave in the k' momentum component.

The integration result can be expressed in terms of the Fourier transform of the density distribution

the term δ((Subscript[E, k'])-(Subscript[E, k])) imposes the energy conservation, which requires |k'|=|k| (the on-shell condition).

Apply to a single atom at the origin: ρ(r)=δ(r),

(OverTilde[ρ])(q)=(∫_r)δ(r)(e^(-iq·r))=(e^(-iq·0))=1.

Therefore the scattering rate is constant for all k' as long as it is on-shell

Γ(k→k')∝δ((Subscript[E, k'])-(Subscript[E, k]))∝δ(|k'|-|k|).

This matches the intuition that an atom scatters a plane wave to a spherical wave.

Fourier Transform of Lattice

(OverTilde[ρ])(q)≡(∫_r)ρ(r)(e^(-iq·r)).

1D lattice (Subscript[R, n])=n a,

Real-space density distribution (modeled by a sequence of Dirac peaks)

ρ(r)=(1/N)(∑_(n∈ℤ))δ(r-(Subscript[R, n]))=(1/N)(∑_(n∈ℤ))δ(r-n a)

Momentum-space density distribution

(OverTilde[ρ])(q)=(1/N)(∑_(n∈ℤ))(e^(-i q (Subscript[R, n])))=(1/N)(∑_(n∈ℤ))(e^(-i q n a))
=((2π)/(N a))(∑_(m∈ℤ))δ(q-((2π m)/a))=((2π)/V)(∑_(m∈ℤ))δ(q-(Subscript[G, m])).

Also a sequence of Dirac peaks located at reciprocal lattice sites:

(Subscript[G, m])=((2π)/a)m.

Argument: constructive interference from all lattice sites happens only when the wave length λ=a/m is an integer fraction of the lattice spacing a.

Plot[Evaluate@Table[Cos[2π r l],{l,0,2}],{r,-2,2},GridLines->{Range[-2,2],None},GridLinesStyle->Directive[Black, AbsoluteThickness[1]],FrameTicks->{{Automatic,Automatic},{Range[-2,2],{#,}&/@Range[-2,2]}},AspectRatio->0.3,FrameLabel->{"\!\(r/a\)",None},ImageSize->{Automatic,120}]

Show that in the functional sense
(i) (∑_(n∈ℤ))(e^(-i n θ))=2π (∑_(m∈ℤ))δ(θ-2π m),
(ii) δ(q a)=(1/a)δ(q).

(∑_(n=-(N-1)/2)^((N-1)/2))(e^(-i n θ))=(e^(i (((N-1))/2) θ))(∑_(n=0)^(N-1))(e^(-i n θ))
=(e^(i ((N θ)/2)))(e^(-((i θ)/2)))((1-(e^(-i N θ)))/(1-(e^(-i θ))))
=(((e^(i ((N θ)/2)))-(e^(-i ((N θ)/2))))/((e^((i θ)/2))-(e^(-((i θ)/2)))))
=((sin ((N θ)/2))/(sin (θ/2))).

(lim_(N→∞))((sin ((N θ)/2))/(sin (θ/2)))(Overscript[=, functional])c (∑_(m∈ℤ))δ(θ-2π m).

Block[{Ν=15},
Plot[Sin[Ν θ/2]/Sin[θ/2],{θ,-4π,4π},PlotStyle->Black,PlotRange->All,PlotLabel->Style[StringTemplate["\!\(N=`1`\)"][Ν],14,Black],FrameLabel->{"θ","\!\(\(sin N θ/2\)\/\(sin θ/2\)\)"},ImageSize->{Automatic,150}]]

(∫_(-π)^π)dθ((sin ((N θ)/2))/(sin (θ/2)))=2π,
(∫_(-π)^π)dθ c(∑_(m∈ℤ))δ(θ-2π m)=(∫_(-π)^π)dθ c δ(θ)=c.

(∑_(n∈ℤ))(e^(-i n θ))=(lim_(N→∞))(∑_(n=-(N-1)/2)^((N-1)/2))(e^(-i n θ))=2π (∑_(m∈ℤ))δ(θ-2π m).

δ(q a)=c δ(q).

∫dq δ(q a)=(1/a)∫d(q a) δ(q a)
=(1/a)∫dx δ(x)
=(1/a),
∫dq c δ(q)=c.

(∑_(n∈ℤ))(e^(-i q n a))=(∑_(n∈ℤ))(e^(-i n (q a)))
=2π (∑_(m∈ℤ))δ(q a-2π m)
=2π (∑_(m∈ℤ))δ((q-((2π m)/a))a)
=((2π)/a)(∑_(m∈ℤ))δ(q-((2π m)/a)).

2D lattices (numerical examples)

Square lattice → square lattice

Block[{L=11,Q=4π,Rs,qs},
Rs=Cases[Tuples[Range[-12,12],2].{{1,0},{0,1}},p_/;Norm[p]<=L];
qs=N@Table[{qx,qy},{qy,Q,-Q,-2π/L},{qx,-Q,Q,2π/L}];
Row@{Graphics[{Black,Point[Rs]},PlotRange->2.5,PlotRangeClipping->True,Frame->True,PlotLabel->Style["\!\(ρ(r)\)",14,Black],ImageSize->{Automatic,140}],ComplexMatrixPlot[Total[Exp[I tTr[{qs,Rs},{{3,5}}]],{-1}]/Length[Rs],DataRange->{{-1,1},{-1,1}}Q,Frame->True,FrameStyle->Black,FrameTicks->{Range[-Q,Q,2π],Range[-Q,Q,2π]},PlotLabel->Style["\!\(ρ\&~(q)\)",14,Black],ImageSize->{Automatic,140}]}]

Block[{L=12,Q=4π,Rs,qs},
Rs=Cases[Tuples[Range[-12,12],2].{{1.5,0},{0,1.5}},p_/;Norm[p]<=L];
qs=N@Table[{qx,qy},{qy,Q,-Q,-2π/L},{qx,-Q,Q,2π/L}];
Row@{Graphics[{Black,Point[Rs]},PlotRange->2.5,PlotRangeClipping->True,Frame->True,PlotLabel->Style["\!\(ρ(r)\)",14,Black],ImageSize->{Automatic,140}],ComplexMatrixPlot[Total[Exp[I tTr[{qs,Rs},{{3,5}}]],{-1}]/Length[Rs],DataRange->{{-1,1},{-1,1}}Q,Frame->True,FrameStyle->Black,FrameTicks->{Range[-Q,Q,2π],Range[-Q,Q,2π]},PlotLabel->Style["\!\(ρ\&~(q)\)",14,Black],ImageSize->{Automatic,140}]}]

Direct lattice expands in the real space ⇔ reciprocal lattice shrinks in the momentum space (hence the name “reciprocal”). Reason: momentum and coordinate scales are reciprocal to each other (q∼(a^(-1))).

Triangular lattice → triangular lattice

Block[{L=12,Q=4π,Rs,rs,qs},
Rs=Cases[Tuples[Range[-14,14],2].{{1,0},{1/2,Sqrt[3]/2}},p_/;Norm[p]<=L];
qs=N@Table[{qx,qy},{qy,Q,-Q,-2π/L},{qx,-Q,Q,2π/L}];
Row@{Graphics[{Point[Rs]},PlotRange->2.2,PlotRangeClipping->True,Frame->True,PlotLabel->Style["\!\(ρ(r)\)",14,Black],ImageSize->{Automatic,140}],ComplexMatrixPlot[Total[Exp[I tTr[{qs,Rs},{{3,5}}]],{-1}]/Length[Rs],DataRange->{{-1,1},{-1,1}}Q,Frame->True,FrameStyle->Black,FrameTicks->{Range[-Q,Q,2π],Range[-Q,Q,2π]},PlotLabel->Style["\!\(ρ\&~(q)\)",14,Black],ImageSize->{Automatic,140}]}]

Reciprocal lattice shares the same point group symmetry as the direct lattice. (In this case, six-fold rotation, horizontal and vertical reflections).

Honeycomb lattice (triangular Bravais lattice + sub-lattice structure) → triangular lattice with non-trivial structure factor.

Block[{L=12,Q=4π,Rs,qs},
Rs=Cases[Join@@Table[(#+r)&/@#,{r,{{Sqrt[3],1},{Sqrt[3],1}/2}}]&[Tuples[Range[-14,14],2].{{Sqrt[3],0},{Sqrt[3]/2,3/2}}],p_/;Norm[p]<=L];
qs=N@Table[{qx,qy},{qy,Q,-Q,-2π/L},{qx,-Q,Q,2π/L}];
Row@{Graphics[{Point[Rs]},PlotRange->2.9,PlotRangeClipping->True,Frame->True,PlotLabel->Style["\!\(ρ(r)\)",14,Black],ImageSize->{Automatic,140}],ComplexMatrixPlot[Total[Exp[I tTr[{qs,Rs},{{3,5}}]],{-1}]/Length[Rs],DataRange->{{-1,1},{-1,1}}Q,Frame->True,FrameStyle->Black,FrameTicks->{Range[-Q,Q,2π],Range[-Q,Q,2π]},PlotLabel->Style["\!\(ρ\&~(q)\)",14,Black],ImageSize->{Automatic,140}]}]

red dot - positive peak, green dot - negative peak. Non-trivial sub-lattice structure leads to non-trivial structure factor (To be discussed in more details).

In general, for atoms at (Subscript[R, n,s]) on a lattice, the scatterer density distribution can be modeled by

where

(Subscript[f, s]) is the scattering amplitude of the type-s atom (different types of atom can scatter the wave with different amplitudes). Normalization: (∑_s)(Subscript[f, s])=1.

N - number of unit cells. Assume thermodynamic limit N->∞. The average  (1/N)(∑_R)… has a well-defined limit as N→∞. The normalization condition is (∫_r)ρ(r)=1.

Therefore by Eq. (),

(OverTilde[ρ])(q)=(∫_r)ρ(r)(e^(-iq·r))
=(1/N)(∫_r)(∑_(n,s))(Subscript[f, s])δ(r-(Subscript[R, n,s]))(e^(-iq·r))
=(1/N)(∑_(n,s))(Subscript[f, s])(e^(-i q·(Subscript[R, n,s]))).

Given (Subscript[R, n,s])=(Subscript[R, n])+(Subscript[r, s]) (recall Eq. ()), the summation can be split into two parts

(OverTilde[ρ])(q)=(1/N)(∑_n)(e^(-i q·(Subscript[R, n])))(∑_s)(Subscript[f, s])(e^(-i q·(Subscript[r, s]))).

If the lattice is a Bravais lattice, there is only one type of atom (no sub-lattice structure), the second summation is just a trivial factor f, which may be rescaled to f=1, such that the scattering amplitude (OverTilde[ρ])(q) is given by the Fourier series

(OverTilde[ρ])(q)=(1/N)(∑_n)(e^(-i q·(Subscript[R, n]))).

Expect (OverTilde[ρ])(q) to peak at multiple points arranged periodically in the momentum space, forming the reciprocal lattice.

Scattering Amplitude

Γ(k→k')∝((|(OverTilde[ρ])(k'-k)|)^2)δ((Subscript[E, k'])-(Subscript[E, k])),

(OverTilde[ρ])(q)=(1/N)(∑_n)(e^(-i q·(Subscript[R, n])))(∑_s)(Subscript[f, s])(e^(-i q·(Subscript[r, s]))).

Given the knowledge of reciprocal lattice, the first summation can be evaluated as (use (Exc. ) results and generalize to higher dimensions)

(∑_n)(e^(-i q·(Subscript[R, n])))=((((2π))^D)/Ω)(∑_m)δ(q-(Subscript[G, m])),

D - space dimension.

Ω - unit cell volume, as defined in Eq. ().

(Subscript[G, m]) - reciprocal lattice vectors, as in Eq. (), labeled by Miller indices m.

If q is a reciprocal lattice point, (e^(-i q·(Subscript[R, n]))) will be unity and the sum is infinite, however if q is not on the reciprocal lattice, the sum will be zero due to the phase cancellation of (e^(-i q·(Subscript[R, n]))). Such behavior is exactly described by the delta function on the right-hand-side.

Then Eq. () can be written as (note N Ω=V is the total volume)

(OverTilde[ρ])(q)=((((2π))^D)/V)(∑_m)δ(q-(Subscript[G, m]))(∑_s)(Subscript[f, s])(e^(-i q·(Subscript[r, s])))
=((((2π))^D)/V)(∑_m)δ(q-(Subscript[G, m]))(∑_s)(Subscript[f, s])(e^(-i (Subscript[G, m])·(Subscript[r, s])))

where we are free to replace q by (Subscript[G, m]) because (OverTilde[ρ])(q) is non-vanishing only when q=(Subscript[G, m]) (when the replacement is valid).

a collection of peaks at q=(Subscript[G, m]) described by δ(q-(Subscript[G, m])),

whose amplitude is given by the structural factor (Subscript[S, m]) (it is determined by the sub-lattice structure, and appears as a factor in front of the delta function peak, hence the name).

Show that Eq. () indeed restores Eq. () given Eq. ().

ρ(r)=(1/(((2π))^D))(∫_q)(OverTilde[ρ])(q)(e^(i q·r))
=(1/V)(∫_q)(∑_m)(Subscript[S, m])δ(q-(Subscript[G, m]))(e^(i q·r))
=(1/V)(∑_m)(Subscript[S, m])(e^(i (Subscript[G, m])·r))
=(1/V)(∑_m)(∑_s)(Subscript[f, s])(e^(-i (Subscript[G, m])·(Subscript[r, s])))(e^(i (Subscript[G, m])·r))
=(1/V)(∑_s)(Subscript[f, s])(∑_m)(e^(i (Subscript[G, m])·(r-(Subscript[r, s])))).

(∑_m)(e^(i (Subscript[G, m])·r))=Ω(∑_n)δ(r-(Subscript[R, n])),

ρ(r)=(Ω/V)(∑_s)(Subscript[f, s])(∑_n)δ(r-(Subscript[r, s])-(Subscript[R, n]))
=(1/N)(∑_(n,s))(Subscript[f, s])δ(r-(Subscript[R, n,s])).

Laue Condition

conservation of crystal momentum (Laue condition)

conservation of energy (for X-ray (Subscript[E, k])=c ℏ|k|, for neutron / electron (Subscript[E, k])=((ℏ^2)/(2m))((|k|)^2))

Manipulate[Block[{b1,b2,Gs,Q=4π},
b1={0,2π};
b2={-2π,0};
Gs=Tuples[Range[-5,5],2].{b1,b2};
Graphics[{{Purple,Circle[{0,0},k]},{Arrowheads[0.06],Red,Arrow[ReIm@{0,k}],Blue,Arrow[ReIm@{0,#}]&/@Cases[k+Exp[I θ]Gs.{1,I},z_/;Abs[Abs[z]-k]<0.2&&Abs[z-k]>0.01]},{AbsolutePointSize[5],Point[ReIm[k+Exp[I θ]Gs.{1,I}]]}},
PlotRange->Q,PlotRangePadding->Scaled[.02],PlotRangeClipping->True,Frame->True,FrameStyle->Black,FrameTicks->{{Range[-Q,Q,2π],None},{Range[-Q,Q,2π],None}},ImageSize->{Automatic,140}]],
{{k,2π},0,4π,ImageSize->Small},{θ,0,π/2,ImageSize->Small},Paneled->False,SaveDefinitions->True]

k - wave number 2π/λ,

θ - rotation angle of the crystal (reciprocal lattice rotates by the same amount as the direct lattice).

Lattice Momentum

Lattice (Bravais lattice) = the intersection of families of lattice planes of all lattice momenta G

ρ(r)=(1/N)(∑_n)δ(r-(Subscript[R, n]))=(1/V)(∑_m)(e^(i (Subscript[G, m])·r)),

see (Exc. ) for derivation.

Take the triangular lattice for example:

Block[{L=2.2,Rs,b1,b2,plt,ms},
Rs=Cases[Tuples[Range[-3,3],2].{{1,0},{1/2,Sqrt[3]/2}},{x_,y_}/;Abs[x]<=2&&Abs[y]<=2];
b1={0,4π/Sqrt[3]};
b2={-2π,2π/Sqrt[3]};
plt[fun_]:=Graphics[{Raster[Table[fun[{x,y}],{y,-L,L,0.1},{x,-L,L,0.1}],{{-L,-L},{L,L}},ColorFunction->ComplexColor],
{AbsoluteThickness[1],Circle[#,Offset@{2,2}]&/@Rs}},BaselinePosition->(Axis->Baseline),ImageSize->70];
ms={{1,0},{0,1},{1,-1}};
Style[Row[{Row[Table[Show[{plt[Cos[m.{b1,b2}.#]&],Graphics[{Arrowheads[0.12],Arrow[{0.1#,1#}&@Normalize[m.{b1,b2}]]}]},PlotLabel->Style[StringTemplate["\!\(\[LeftBracketingBar]G\_(`1`)⟩\)"][Row[m,","]],14,Black]],{m,ms}],"+"],"=",plt[Sum[Cos[m.{b1,b2}.#],{m,ms}]/Sqrt[3]&]}],14,FontFamily->"LaTeX"]
]

Each family of lattice planes is in one-to-one correspondence with the possible directions of reciprocal lattice vectors G, to which they are normal.

The length of G is always an integer multiple of 2π/d,

|G|=((2π m)/d),

where d is the spacing between the lattice planes. (The corresponding wave length λ=d/m).

The particle that carries the momentum ℏ G is the condensed phonon in the crystal (To be explained later).

Bragg Condition

An incoming wave is reflected off of a family of adjacent layers of atoms separated by the inter-layer distance d.

Block[{θ=π/3,R=2.5},Row[{Graphics[{Table[{{AbsoluteThickness[1],Dashing[0.04],Gray,Line[{{-2,y},{2,y}}]},{Arrowheads[{{0.01,0.8,Graphics[Line[Offset/@{{-5,-3},{0,0},{-5,3}}]]}}],Arrow[{{-R Cos[θ],R Sin[θ]}+y Cos[θ]{Sin[θ],Cos[θ]},{0,y}}],Arrow[{{0,y},{R Cos[θ],R Sin[θ]}+y Cos[θ]{-Sin[θ],Cos[θ]}}]},{AbsolutePointSize[5],Point[Table[{x,y},{x,-2,2}]]}},{y,0,-2,-1}],
Translate[{AbsoluteThickness[1],Arrowheads[{-0.06,0.06}],Arrow[{0,#}&/@{0,-1}],Text["\!\( d\)",{0,-0.5},{-1,0}]},{-1.8,-1}],
Translate[{AbsoluteThickness[1],{Dashing[{0.02,0.02}],Line[{{0,0},2Cos[θ]{Sin[θ],-Cos[θ]}}],Line[{{0,-1},Cos[θ]{Sin[θ],-Cos[θ]}+{0,-1}}]},
Translate[{Rotate[Line[Offset/@{{-4,0},{-4,4},{0,4}}],θ,{0,0}]},Cos[θ]{Sin[θ],-Cos[θ]}],Translate[{Arrowheads[{-0.06,0.06}],Arrow[{{0,0},Sin[θ]{Cos[θ],Sin[θ]}}],Text["\!\(d sin θ\)",Sin[θ]/2{Cos[θ],Sin[θ]},{-1.2,0}]},0.5Cos[θ]{Sin[θ],-Cos[θ]}+{0,-1}]},{0,-1}],
Translate[{AbsoluteThickness[1],Circle[{0,0},0.3,{π-θ,π}],Text["θ",ReIm[0.3 Exp[I (π-θ/2)]],{1,0}]},{0,-2}]},ImageSize->120],
Graphics[{{Arrowheads[0.4],Arrow[{{0,0},{Cos[θ],-Sin[θ]}}],Arrow[{{0,0},{Cos[θ],Sin[θ]}}],Arrow[{{Cos[θ],-Sin[θ]},{Cos[θ],Sin[θ]}}]},
{AbsoluteThickness[1],{Dashed,Line[{{0,0},{Cos[θ],0}}]},Circle[{0,0},0.2,{0,θ}],Text["θ",{0.3,0.15}]},
{Text["\!\( k\)",{Cos[θ],-Sin[θ]}/2,{1,1}],
Text["\!\( k\^′\)",{Cos[θ],Sin[θ]}/2,{0.2,-1}],
Text["\!\( G\)",{Cos[θ],0},{-1,0}]}},ImageSize->36]},Spacer[10]]]

The wave has been deflected by angle 2θ.

Wave deflected from the next deeper layer travels additional distance 2d sin θ (extra d sin θ accumulated on both the incoming and outgoing rays), compare to that from the previous layer.

Condition for constructive interference (Bragg condition)

with m∈ℤ.

2 d sin θ=m λ⇔k'-k=G.

Direction: the direction of k'-k is always normal to the lattice plane, matching the direction of G.

Length: |k'|=|k|=2π/λ, |G|=2π m/d.

2((2π)/λ)sin θ=((2π m)/d)=>2d sin θ=m λ.

Structural Factor

(Subscript[S, m])=(∑_s)(Subscript[f, s])(e^(-i (Subscript[G, m])·(Subscript[r, s]))),

(Subscript[f, s]) is the atomic form factor (scattering amplitude) of the type-s atom,

(Subscript[r, s]) is the sub-lattice vector of the type-s atom.

(Subscript[S, m])=(∑_s)(Subscript[f, s])exp(-i (∑_(i=1)^D)(Subscript[m, i])(Subscript[b, i])·(∑_(j=1)^D)(Subscript[n, s,j])(Subscript[a, j]))
=(∑_s)(Subscript[f, s])exp(-i (∑_(i,j))(Subscript[m, i])(Subscript[n, s,j])(Subscript[b, i])·(Subscript[a, j]))
=(∑_s)(Subscript[f, s])exp(-2π i (∑_(i,j))(Subscript[m, i])(Subscript[n, s,j])(KroneckerDeltaSeq[i,j]))
=(∑_s)(Subscript[f, s])exp(-2π i (∑_i)(Subscript[m, i])(Subscript[n, s,i])),

Honeycomb lattice (recall Eq. ())

({{{, ({{(Subscript[a, 1])=((Sqrt[3]),0)}, {(Subscript[a, 2])=(((Sqrt[3])/2),(3/2))}})}}),({{{, ({{(Subscript[r, A])=((Sqrt[3]),1)}, {(Subscript[r, B])=(((Sqrt[3])/2),(1/2))}})}}).

Block[{L=1,a1,a2,Rs},
a1={Sqrt[3],0};
a2={Sqrt[3],3}/2;
Rs=Tuples[Range[-L-2,L+2],2].{a1,a2};
Graphics[{{LightGray,AbsoluteThickness[1],Translate[Line[{{0,0},#}]&/@CirclePoints[{1,π/2},3],#+{Sqrt[3],1}]&/@Rs},{AbsolutePointSize[5],Blue,Point[#+{Sqrt[3]/2,1/2}&/@Rs],Red,Point[#+{Sqrt[3],1}&/@Rs]},{Arrowheads[0.06],{Arrow[{{0,0},#1}],Text[StringTemplate["\!\(\*StyleBox[\"a\",Bold]\_`1`\)"][#2],#1,#3]}&@@@{{a1,"1",{0,0.8}},{a2,"2",{-1.,-0.2}}}},
{Blue,Text["\!\(,B\)",{Sqrt[3]/2,1/2},{1.5,0.5}],Red,Text["\!\(,A\)",{Sqrt[3],1},{-1.5,-0.5}]}},PlotRange->2L,PlotRangeClipping->True,PlotRangePadding->0.8,Frame->True,FrameTicks->{{Range[-2L,2L],{#,""}&/@Range[-2L,2L]},{Range[-2L,2L],{#,""}&/@Range[-2L,2L]}},ImagePadding->{{Automatic, 5}, {Automatic, 5}},ImageSize->140]]

Rewrite (Subscript[r, s]) to (Subscript[n, s]) by Eq. (),

({{{, ({{(Subscript[n, A])=((2/3),(2/3))}, {(Subscript[n, B])=((1/3),(1/3))}})}}).

The structure factor is given by Eq. ()

(Subscript[S, m])=(Subscript[f, A])(e^(-2π i 2((Subscript[m, 1])+(Subscript[m, 2]))/3))+(Subscript[f, B])(e^(-2π i ((Subscript[m, 1])+(Subscript[m, 2]))/3))
=(Subscript[f, A])(ϖ^((Subscript[m, 1])+(Subscript[m, 2])))+(Subscript[f, B])(ϖ^(-((Subscript[m, 1])+(Subscript[m, 2])))).

where ϖ=(e^(2π i/3 )) is the third root of unity.

If (Subscript[f, A])=(Subscript[f, B])=f (A, B sub-lattices are equivalent as in graphene)

(Subscript[S, m])=({{{, ({{2f, (Subscript[m, 1])+(Subscript[m, 2])(Overscript[=, mod 3])0}, {-f, (Subscript[m, 1])+(Subscript[m, 2])(Overscript[=, mod 3])1,2}})}}).

The reciprocal lattice is given by (Subscript[G, m])=(Subscript[m, 1])(Subscript[b, 1])+(Subscript[m, 2])(Subscript[b, 2]) with

({{{, ({{(Subscript[b, 1])=(((2π)/(Sqrt[3])),-((2π)/3))}, {(Subscript[b, 2])=(0,((4π)/3))}})}}).

Block[{b1,b2,Q=4π},
b1={2π/Sqrt[3],-2π/3};
b2={0,4π/3};
Graphics[{{Arrowheads[0.1],Arrow[{{0,0},#}]&/@{b1,b2}},{Text["\!\(b\_"<>#1<>"\)",##2]&@@@{{"1",b1,{-1.2,0.3}},{"2",b2,{0,-0.7}}}},{If[Mod[#1+#2,3]==0,Red,Lighter[Green]],Point[#1 b1+#2 b2]}&@@@Tuples[Range[-4,4],2]},PlotRange->Q,PlotRangePadding->Scaled[.02],PlotRangeClipping->True,Frame->True,FrameStyle->Black,FrameTicks->{{Range[-Q,Q,2π],None},{Range[-Q,Q,2π],None}},ImageSize->{Automatic,120}]]

red - (Subscript[m, 1])+(Subscript[m, 2])(Overscript[=, mod 3])0 => (Subscript[S, m])=2f,

light green - (Subscript[m, 1])+(Subscript[m, 2])(Overscript[=, mod 3])1,2 => (Subscript[S, m])=-f.

Matching the numerical result of lattice Fourier transform.

Block[{L=12,Q=4π,Rs,qs},
Rs=Cases[Join@@Table[(#+r)&/@#,{r,{{Sqrt[3],1},{Sqrt[3],1}/2}}]&[Tuples[Range[-14,14],2].{{Sqrt[3],0},{Sqrt[3]/2,3/2}}],p_/;Norm[p]<=L];
qs=N@Table[{qx,qy},{qy,Q,-Q,-2π/L},{qx,-Q,Q,2π/L}];
ComplexMatrixPlot[Total[Exp[I tTr[{qs,Rs},{{3,5}}]],{-1}]/Length[Rs],DataRange->{{-1,1},{-1,1}}Q,Frame->True,FrameStyle->Black,FrameTicks->{Range[-Q,Q,2π],Range[-Q,Q,2π]},PlotLabel->Style["\!\(ρ\&~(q)\)",14,Black],ImageSize->{Automatic,140}]]

Cubic lattice (recall Eq. ()) with both corner (A) and center (B) sites

({{{, ({{(Subscript[a, 1])=(1,0,0)}, {(Subscript[a, 2])=(0,1,0)}, {(Subscript[a, 3])=(0,0,1)}})}}), ({{{, ({{(Subscript[r, A])=(0,0,0)}, {(Subscript[r, B])=((1/2),(1/2),(1/2))}})}}).

Block[{pts},
pts={{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1},{1,1,1}/2};
Graphics3D[{GraphicsComplex[pts,
{{Red,Sphere[#,0.08]&/@Range[8],Text["\!\( A\)",2,{-1,1}]},{Blue,Sphere[9,0.08],Text["\!\( B\)",9,{0,-1}]},{Black,Tube[{{1,2,4,3,1},{5,6,8,7,5},{1,5},{2,6},{3,7},{4,8}}]},{LightGray,Tube[{9,#}&/@Range[8]]},{Black,Arrowheads[{{0.1,0.95}}],Arrow[Tube[{1,#},0.02]]&/@{2,3,5}},{Text[StringTemplate["\!\(\*StyleBox[\"a\",Bold]\_`1`\)"][#2],#1,#3]&@@@{{2,1,{1,-1}},{3,2,{1,-1}},{5,3,{-1.5,1}}}}}]},Lighting->"Neutral",ViewPoint->{2,0,1},Boxed->False,ImageSize->120]]

Rewrite (Subscript[r, s]) to (Subscript[n, s]) by Eq. (),

({{{, ({{(Subscript[n, A])=(0,0,0)}, {(Subscript[n, B])=((1/2),(1/2),(1/2))}})}}).

The structure factor is given by Eq. ()

(Subscript[S, m])=(Subscript[f, A])+(Subscript[f, B])(e^(-2π i ((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3]))/2))
=(Subscript[f, A])+(((Subscript[f, B])(-1))^((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3]))).

If (Subscript[f, A])=(Subscript[f, B])=f (A, B sub-lattices are equivalent as in body-centered cubic (BCC) lattice)

(Subscript[S, m])=({{{, ({{2f, (Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3]) is even}, {0, (Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3]) is odd}})}}).

The reciprocal lattice is given by (Subscript[G, m])=(∑_(i=1)^3)(Subscript[m, i])(Subscript[b, i]) with

({{{, ({{(Subscript[b, 1])=2π (1,0,0)}, {(Subscript[b, 2])=2π (0,1,0)}, {(Subscript[b, 3])=2π (0,0,1)}})}}).

Block[{b1,b2,b3},
{b1,b2,b3}=2π IdentityMatrix[3];
Graphics3D[{Table[{If[EvenQ[m1+m2+m3],Black,Opacity[0.5]],Sphere[m1 b1+m2 b2+m3 b3,0.5]},{m1,0,2},{m2,0,2},{m3,0,2}],
{Black,Rotate[Translate[Translate[Tube[{{0,0,0},2b1},0.05],# b2]&/@Range[0,2],# b3]&/@Range[0,2],2π #/3,{1,1,1},{0,0,0}]&/@Range[3]}},Lighting->"Neutral",ViewPoint->{3,-1,1},Boxed->False,ImageSize->140]]

Solid (black) - (Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3]) even => (Subscript[S, m])≠0. They form a face-centered cubic (FCC) lattice, verifying the fact that the reciprocal lattice of BCC lattice is a FCC lattice (see (HW )).

Transparent - (Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3]) odd => (Subscript[S, m])=0. These are called extinction points in the diffraction pattern, where the scattering wave is systematically absent.

Symmetry reason for the extinction phenomenon

When A, B sites are equivalent, the lattice symmetry is enlarged. For example, translation by the vector (Subscript[r, B])=((1/2),(1/2),(1/2)) becomes a symmetry (which was not for a simple cubic lattice). This means

ρ(r)=ρ(r+(Subscript[r, B])),

is a new property of the density distribution ρ.

Under the symmetry action, the scattering amplitude (OverTilde[ρ])(q) transforms as (according to Eq. ())

(OverTilde[ρ])(q)=(∫_r)ρ(r)(e^(-iq·r))
→(∫_r)ρ(r+(Subscript[r, B]))(e^(-iq·(r+(Subscript[r, B]))))
=(e^(-iq·(Subscript[r, B])))(∫_r)ρ(r)(e^(-iq·r))
=(e^(-iq·(Subscript[r, B])))(OverTilde[ρ])(q).

Apply to the reciprocal lattice points

(OverTilde[ρ])((Subscript[G, m]))→(e^(-(Subscript[iG, m])·(Subscript[r, B])))(OverTilde[ρ])((Subscript[G, m]))
=(e^(-2π i m·(Subscript[n, B])))(OverTilde[ρ])((Subscript[G, m]))
=(((-1))^((Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3])))(OverTilde[ρ])((Subscript[G, m])).

Therefore, if (Subscript[m, 1])+(Subscript[m, 2])+(Subscript[m, 3]) is odd, Eq. () implies that (OverTilde[ρ])((Subscript[G, m])) will transform to -(OverTilde[ρ])((Subscript[G, m])). The only way to preserve (OverTilde[ρ])((Subscript[G, m])) under the symmetry transformation is to set (OverTilde[ρ])((Subscript[G, m]))=0, hence the extinction phenomenon.

Summary

Relation between lattice and reciprocal lattice.

({{Bravais lattice,  , Reciprocal lattice}, {(Subscript[R, n])=(∑_(i=1)^D)(Subscript[n, i])(Subscript[a, i]), (Underscript[⇋, (Subscript[b, i])·(Subscript[a, j])=2π (Subscript[δ, i,j])]), (Subscript[G, m])=(∑_(i=1)^D)(Subscript[m, i])(Subscript[b, i])}, {(↓^( sub-lattice)),  ,  }, {(Subscript[R, n,s])=(Subscript[R, n])+(Subscript[r, s]),  ,  }, {General lattice,  ,  }})

Volume of primitive unit cell

Ω=|(Subscript[a, 1])·((Subscript[a, 2])×(Subscript[a, 3]))|.

Construct (Subscript[b, i]) from (Subscript[a, i])

(Subscript[b, i])=((2π)/Ω)(1/2)(Subscript[ϵ, i,j,k]) (Subscript[a, j])×(Subscript[a, k]).

Fourier transformation between atomic density distribution ρ(r) and scattering amplitude (OverTilde[ρ])(q)

({{ρ(r)=(1/N)(∑_(n,s))(Subscript[f, s])δ(r-(Subscript[R, n,s])), (Underoverscript[⇌, (ρ(r)=(1/(((2π))^D))(∫_q)(OverTilde[ρ])(q)(e^(i q·r))), ((OverTilde[ρ])(q)=(∫_r)ρ(r)(e^(-iq·r)))]), (OverTilde[ρ])(q)=((((2π))^D)/V)(∑_m)(Subscript[S, m])δ(q-(Subscript[G, m]))}, {  ,  , ((Subscript[S, m])=(∑_s)(Subscript[f, s])(e^(-i (Subscript[G, m])·(Subscript[r, s]))))}})

(OverTilde[ρ])(q) determines the scattering rate

Γ(k→k')∝((|(OverTilde[ρ])(k'-k)|)^2)δ((Subscript[E, k'])-(Subscript[E, k])).

Diffraction spot location is determined by either the Laue condition or the Bragg condition (they are equivalent)

({{{, ({{k'-k=(Subscript[G, m])}, {|k'|=|k|}})}})⇔2 d sin θ=m λ.

Diffraction spot intensity ((|(Subscript[S, m])|)^2)  is determined by the structural factor

(Subscript[S, m])=(∑_s)(Subscript[f, s])(e^(-2πi m·(Subscript[n, s]))).

Diffraction spot can extinct when (Subscript[S, m])=0.

Phonon

Lattice Vibration

Hamiltonian of Atoms

(Subscript[p, i]) - the momentum of the ith atom,

(Subscript[r, i]) - the position of the ith atom,

(Subscript[m, i]) - the mass of the ith atom.

V((Subscript[r, 1]),(Subscript[r, 2]),…)=(∑_(i≠j))V((Subscript[r, i])-(Subscript[r, j])),

Show[Plot[1/r^12-2/r^6,{r,0.8,1.6},PlotStyle->Black,PlotRange->{-1.2,1.},FrameLabel->{"\!\(r/a\)","\!\(V(r)/V\_0\)"},GridLines->{{1},None},ImageSize->200],Plot[-1+36(r-1)^2,{r,0.9,1.1},PlotStyle->Red]]

V(r)=(Subscript[V, 0])+(κ/2)(((r-a))^2)+…,

Harmonic Approximation

Let (Subscript[R, i]) be the equilibrium position of the ith atom, and (Subscript[x, i]) be the displacement of the ith atom from its equilibrium position, s.t.

(Subscript[r, i])=(Subscript[R, i])+(Subscript[x, i]).

Assuming (Subscript[x, i]) are small, expand the potential energy to the 2nd order of (Subscript[x, i])

V((Subscript[r, 1]),(Subscript[r, 2]),…)=V((Subscript[R, 1])+(Subscript[x, 1]),(Subscript[R, 2])+(Subscript[x, 2]),…)
=V((Subscript[R, 1]),(Subscript[R, 2]),…)
+(∑_i)(Subscript[x, i])·(Subscript[∇, (Subscript[R, i])])V((Subscript[R, 1]),(Subscript[R, 2]),…)
+(1/2)(∑_(i,j))((Subscript[x, i])·(Subscript[∇, (Subscript[R, i])]))((Subscript[x, j])·(Subscript[∇, (Subscript[R, j])]))V((Subscript[R, 1]),(Subscript[R, 2]),…)
+…

where

The 0th order term V((Subscript[R, 1]),(Subscript[R, 2]),…)=(Subscript[V, eq]) is a constant shift of energy, which can be dropped from the Hamiltonian without affecting dynamics.

The 1st order term vanishes, because the gradient (Subscript[∇, (Subscript[R, i])])V=0 vanishes at equilibrium (by definition).

The 2nd order term takes the form of

(1/2)(∑_(i,j))(Subscript[x, i])·(Subscript[K, i,j])·(Subscript[x, j]),

where (Subscript[K, i,j]) is a matrix for every i,j, and is called the elastic tensor,

(Subscript[K, i,j])=(Subscript[∇, (Subscript[R, i])])(Subscript[∇, (Subscript[R, j])])V((Subscript[R, 1]),(Subscript[R, 2]),…).

Apply to the pairwise interaction model Eq. (), and assume the inter-atomic interaction V(r)=V(r) is isotropic, by Taylor expansion,

(1/2)(∑_(i,j))(Subscript[x, i])·(Subscript[K, i,j])·(Subscript[x, j])=(1/2)(∑_(i≠j))((((Subscript[κ, i,j]^⊥)((Subscript[x, i])-(Subscript[x, j])))^2)+((Subscript[κ, i,j]^∥)-(Subscript[κ, i,j]^⊥))((((Subscript[(OverHat[R]), i,j])·((Subscript[x, i])-(Subscript[x, j]))))^2) ),

where

(Subscript[(OverHat[R]), i,j])=(Subscript[R, i,j])/|(Subscript[R, i,j])| is the unit vector along the bond displacement (Subscript[R, i,j])=(Subscript[R, i])-(Subscript[R, j]),

Graphics[{Arrowheads[0.1],Arrow[{{0,0},{4,3}}],Arrow[{{0,0},{5,0}}],{Dashed,Line[{{4,0},{4,3}}]},Text["\!\(δx\_\(i,j\)=x\_i-x\_j\)",{2,1.5},{0,-1},{4,3}],
Text["\!\(R\&^\_\(i,j\)·δx\_\(i,j\)\)",{2,0},{0,1}],
Text["\!\(R\&^\_\(i,j\)\)",{5,0},{-1,0}],
Text["\!\(\@\(δx\_\(i,j\)\%2-(R\&^\_\(i,j\)·δ\!x\_\(i,j\))\^2\)\)",{4,1.5},{-1.1,0}]},ImageSize->{Automatic,95}]

(Subscript[κ, i,j]^∥) and (Subscript[κ, i,j]^⊥) are elastic constants along the parallel and perpendicular directions of the bond, which are defined by

(Subscript[κ, i,j]^∥)=V''((Subscript[R, i,j])), (Subscript[κ, i,j]^⊥)=((V'((Subscript[R, i,j])))/(Subscript[R, i,j])).

Derive Eq. ().

V((Subscript[r, 1]),(Subscript[r, 2]),…)=(∑_(i≠j))V((Subscript[r, i])-(Subscript[r, j])),

V((Subscript[r, 1]),(Subscript[r, 2]),…)=V((Subscript[R, 1])+(Subscript[x, 1]),(Subscript[R, 2])+(Subscript[x, 2]),…)
=(∑_(i≠j))V((Subscript[R, i])-(Subscript[R, j])+(Subscript[x, i])-(Subscript[x, j]))
=(∑_(i≠j))V((Subscript[R, i,j])+(Subscript[x, i,j])),

|(Subscript[R, i,j])+(Subscript[x, i,j])|=(Sqrt[(Subscript[R, i,j])·(Subscript[R, i,j])+2(Subscript[R, i,j])·(Subscript[x, i,j])+(Subscript[x, i,j])·(Subscript[x, i,j])])

V((Subscript[r, 1]),(Subscript[r, 2]),…)=(∑_(i≠j))V((Sqrt[(Subscript[R, i,j])·(Subscript[R, i,j])+2(Subscript[R, i,j])·(Subscript[x, i,j])+(Subscript[x, i,j])·(Subscript[x, i,j])])).

V((Subscript[r, 1]),(Subscript[r, 2]),…)=(∑_(i≠j))V((Sqrt[(Subscript[R, i,j])·(Subscript[R, i,j])+2λ (Subscript[R, i,j])·(Subscript[x, i,j])+(λ^2)(Subscript[x, i,j])·(Subscript[x, i,j])]))
=…+((λ^2)/2)(∑_(i≠j))((((((Subscript[R, i,j])·(Subscript[R, i,j]))((Subscript[x, i,j])·(Subscript[x, i,j]))-((((Subscript[R, i,j])·(Subscript[x, i,j])))^2))V'((Subscript[R, i,j])))/(Subscript[R, i,j]^3))+((((((Subscript[R, i,j])·(Subscript[x, i,j])))^2)V''((Subscript[R, i,j])))/(Subscript[R, i,j]^2)))+…,

V((Subscript[r, 1]),(Subscript[r, 2]),…)=(1/2)(∑_(i≠j))((((Subscript[x, i,j])·(Subscript[x, i,j]))-((((Subscript[(OverHat[R]), i,j])·(Subscript[x, i,j])))^2))(((V'((Subscript[R, i,j])))/(Subscript[R, i,j])))+((((Subscript[(OverHat[R]), i,j])·(Subscript[x, i,j])))^2)V''((Subscript[R, i,j])))+….

(Subscript[κ, i,j]^∥)=V''((Subscript[R, i,j])), (Subscript[κ, i,j]^⊥)=((V'((Subscript[R, i,j])))/(Subscript[R, i,j])),

V((Subscript[r, 1]),(Subscript[r, 2]),…)=(1/2)(∑_(i≠j))(((Subscript[x, i,j])·(Subscript[x, i,j]))(Subscript[κ, i,j]^⊥)+((((Subscript[(OverHat[R]), i,j])·(Subscript[x, i,j])))^2) ((Subscript[κ, i,j]^∥)-(Subscript[κ, i,j]^⊥)))+….

Phonon Dispersion Relation

Block[{L=8,r=0.1,x},
x[i_]:=-0.5Cos[i π/6];
Graphics[{{LightGray,Line[{#,0}&/@{0.2,L+0.8}]},Translate[{{{EdgeForm[LightGray],White,Disk[{0,0},r]},{EdgeForm[Black],Disk[{x[#],0},r]}},{AbsoluteThickness[1],Line[{{0,-1.5r},{0,-2.5r}}],If[x[#]≠0,Line[{{{0,-2r},{x[#],-2r}},{{x[#]-0.5r Sign[x[#]],-1.5r},{x[#],-2r},{x[#]-0.5r Sign[x[#]],-2.5r}}}]]}},{#,0}]&/@Range[L],
Text["\!\(x\_\("<>#1<>"\)\)",{#2+x[#2]/2,-0.5}]&@@@{{"i-1",5},{"i",6},{"i+1",7}},
Translate[{AbsoluteThickness[1],Line[{{{0,0},{1,0}},{{0,-0.5r},{0,0.5r}},{{1,-0.5r},{1,0.5r}}}],Text["\!\(,a\)",{0.5,0},{0,-0.5}]},{3,2.5r}]},ImageSize->300]]

White dots - equilibrium positions (Subscript[R, i])=i a.

Black dots - actual positions of atoms (the ith atom displaced by (Subscript[x, i])).

H=(∑_i)((Subscript[p, i]^2)/(2m))+(∑_i)(κ/2)((((Subscript[x, i])-(Subscript[x, i+1])))^2).

(Subscript[K, i,j])=κ(2(KroneckerDeltaSeq[i,j])-(KroneckerDeltaSeq[i,j+1])-(KroneckerDeltaSeq[i,j-1])).

Equation of motion:

m(Subscript[(Overscript[x, ..]), i])=-(∑_j)(Subscript[K, i,j])(Subscript[x, j]).

Derive the equation of motion Eq. () from the Hamiltonian Eq. ().

(Subscript[(OverDot[x]), i])=((∂H)/(∂(Subscript[p, i])))=((Subscript[p, i])/m),
 (Subscript[(OverDot[p]), i])=-((∂H)/(∂(Subscript[x, i])))=-(∑_j)(Subscript[K, i,j])(Subscript[x, j]),

(Subscript[(Overscript[x, ..]), i])=((Subscript[(OverDot[p]), i])/m)=-(1/m)(∑_j)(Subscript[K, i,j])(Subscript[x, j]).

Take the trial solution

(Subscript[x, i])=u (e^(-i ω t+i q (Subscript[R, i]))),

substitute into Eq. (),

-m (ω^2)u (e^(-i ω t+i q (Subscript[R, i])))=-(∑_j)(Subscript[K, i,j])u (e^(-i ω t+i q (Subscript[R, j])))

or

m (ω^2)=(∑_j)(Subscript[K, i,j])(e^(-i q ((Subscript[R, i])-(Subscript[R, j]))))
=κ(∑_j)(2(KroneckerDeltaSeq[i,j])-(KroneckerDeltaSeq[i,j+1])-(KroneckerDeltaSeq[i,j-1]))(e^(-i q a (i-j)))
=κ (2-(e^(-i q a ))-(e^(+i q a )))
=2κ (1-cos q a)
=κ (((2 sin ((q a)/2)))^2).

Plot[2 Abs[Sin[q/2]],{q,-π,π},PlotStyle->Black,Epilog->{{AbsoluteThickness[1],Dashed,Line[{{0,0},{2,2}}],Text["\!\(v\_0 q\)",{1.2,1.2},{0,-0.8},{1,2}]}},FrameLabel->{"\!\( q\)","\!\(ω\_q\)"},FrameTicks->{{Automatic,Automatic},{{π#1,If[Mod[#1,1]==0,Round[#1]"\!\(π/a\)",""],#3}&@@@Charting`ScaledTicks[{Identity,Identity}][-1,1],Automatic}},ImageSize->180]

Each point on the dispersion relation corresponds to a sound wave mode (a normal mode of the 1D lattice vibration). The sound wave quantum is called a phonon. The dispersion relation indicates that a phonon excitation of momentum ℏ q will carry energy ℏ (Subscript[ω, q]).

It is sufficient to describe the dispersion relation within the first Brillouin zone [-π/a,π/a).

Because (Subscript[ω, q])=(Subscript[ω, q±(2π/a)]) is periodic in the momentum space (with the same periodicity as the reciprocal lattice), the behavior of (Subscript[ω, q]) within the 1st BZ fully determines its behavior outside the 1st BZ.

Show[{Plot[2 Abs[Sin[q/2]],{q,-3π,3π},PlotStyle->Directive[Dashed,Gray]],Plot[2 Abs[Sin[q/2]],{q,-π,π},PlotStyle->Black]},FrameLabel->{"\!\( q\)","\!\(ω\_q\)"},PlotRangePadding->{0,0.1},PlotRangeClipping->False,ImagePadding->{{Automatic,Automatic},{Automatic,20}},Prolog->{{Lighter[#1,0.85],Rectangle[{#2,-0.1},{#3,3}]}&@@@Thread[{{Green,Red,Blue,Blue,Red,Green},Most[#],Rest[#]}&@(π Range[-3,3])]},Epilog->{Translate[{#1 ,Text[#2,{0,0},{0,-0.6}]},{#3,2.2}]&@@@{{Blue,"1st BZ",0},{Red,"2nd",-3π/2},{Red,"2nd",3π/2},{Green,"3rd",-5π/2},{Green,"3rd",5π/2}}},FrameTicks->{{Automatic,Automatic},{{π#1,If[Mod[#1,1]==0,Round[#1]"\!\(π/a\)",""],#3}&@@@Charting`ScaledTicks[{Identity,Identity}][-3,3],Automatic}},AspectRatio->1/4,ImageSize->300]

Because q±(2π m/a) describe the same wave mode as q on the lattice (known as the aliasing of waves), there is no need to go beyond the 1st BZ.

Block[{q=π/Sqrt[5],ϕ=π/Sqrt[10]},
Show[{Plot[{Cos[q i+ϕ],Cos[(q+2π)i+ϕ]},{i,-0.5,10.5},PlotStyle->(Lighter[#,0.5]&/@{Blue,Red})],DiscretePlot[Cos[q i+ϕ],{i,0,10},PlotStyle->Directive[Black,AbsolutePointSize[5]]]},AspectRatio->1/4,Axes->False,GridLines->{Range[0,10],{0}},FrameLabel->{"\!\( i\)","\!\(x\_i\)"},FrameTicks->{{Automatic,Automatic},{Range[0,10],None}},ImageSize->300]]

Low-frequency limit (BZ center): linear dispersion (Subscript[ω, q])=(Subscript[v, 0])|q| => constant speed of sound

High-frequency limit (BZ boundary):

Group velocity vanishes at BZ boundary q=±π/a.

Inversion symmetry requires (Subscript[ω, q])=(Subscript[ω, -q]) to be even in q, so its first order derivative (Subscript[v, q])=(∂_q)(Subscript[ω, q]) must be odd in q, i.e. (Subscript[v, q])=-(Subscript[v, -q]). Further consider the periodicity (Subscript[v, q])=(Subscript[v, q-(2π/a)]), we have (Subscript[v, q-(2π/a)])=-(Subscript[v, -q]). Set q=π/a at the BZ boundary, this implies (Subscript[v, -(π/a)])=-(Subscript[v, -(π/a)]), thus (Subscript[v, -(π/a)])=0 (if (Subscript[v, q]) is smooth at the BZ boundary). Similar argument applies to (Subscript[v, (π/a)])=0.

Bragg reflection: q=±π/a are related by a lattice momentum G=(2π)/a, strongly scattered to each other by the lattice => standing wave mode (equal-weight superposition of left- and right-moving modes):

(e^(i (π/a)(Subscript[R, i])))+(e^(-i (π/a)(Subscript[R, i])))∝(((-))^i).

Acoustic and Optical Modes

Block[{L=8,r=0.1,x},
x[i_]:=-0.5Cos[i π/6];
Graphics[{{LightGray,Line[{#,0}&/@{0.2,L+0.8}]},Translate[{{{EdgeForm[LightGray],White,Disk[{0,0},r(1+0.2(-1)^#)]},{EdgeForm[Black],Disk[{x[#],0},r(1+0.2(-1)^#)]}},{AbsoluteThickness[1],Line[{{0,-1.5r},{0,-2.5r}}],If[x[#]≠0,Line[{{{0,-2r},{x[#],-2r}},{{x[#]-0.5r Sign[x[#]],-1.5r},{x[#],-2r},{x[#]-0.5r Sign[x[#]],-2.5r}}}]]}},{#,0}]&/@Range[L],
Text["\!\(x\_\("<>#1<>"\)\)",{#2+x[#2]/2,-0.5}]&@@@{{"i,A",5},{"i,B",6}},
Translate[{AbsoluteThickness[1],Line[{{{0,0},{1,0}},{{0,-0.5r},{0,0.5r}},{{1,-0.5r},{1,0.5r}}}],Text["\!\(,a=1\)",{0.5,0},{0,-0.5}]},{3,2.5r}]},ImageSize->300]]

(Subscript[R, i,A])=2i,
(Subscript[R, i,B])=2i+1.

H=(∑_i)(((Subscript[p, i,A]^2)/(2(Subscript[m, A])))+((Subscript[p, i,B]^2)/(2(Subscript[m, B]))))+(κ/2)(∑_i)(((((Subscript[x, i,A])-(Subscript[x, i,B])))^2)+((((Subscript[x, i,B])-(Subscript[x, (i+1),A])))^2)).

Introduce the following vectors in each unit cell (taking A,B as basis)

The Hamiltonian Eq. () can be generally written as

where M and K matrices are

(Subscript[((M^(-1))), i,j])=(({{(Subscript[m, A]^(-1))(Subscript[δ, i,j]), 0}, {0, (Subscript[m, B]^(-1))(Subscript[δ, i,j])}})),
(Subscript[K, i,j])=κ(({{2(Subscript[δ, i,j]), -(Subscript[δ, i,j])-(Subscript[δ, i,j+1])}, {-(Subscript[δ, i,j])-(Subscript[δ, i,j-1]), 2(Subscript[δ, i,j])}})).

(∑_j)(Subscript[M, i,j])(Subscript[(Overscript[x, ..]), j])=-(∑_j)(Subscript[K, i,j])(Subscript[x, j]).

Take the trial solution (labeled by q)

(Subscript[x, i])=(({{(Subscript[x, i,A])}, {(Subscript[x, i,B])}}))=(e^(-i ω t))(({{(Subscript[u, A])(e^(i q (Subscript[R, i,A])))}, {(Subscript[u, B])(e^(i q (Subscript[R, i,B])))}})),

substitute into Eq. (),

(ω^2)(({{(Subscript[m, A]), 0}, {0, (Subscript[m, B])}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}}))=2κ(({{1, -cos q}, {-cos q, 1}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}}))

Derive Eq. ().

lhs=(∑_j)(Subscript[M, i,j])(Subscript[(Overscript[x, ..]), j])=(∑_j)(Subscript[M, i,j])(∂_t^2)(e^(-i ω t))(({{(Subscript[u, A])(e^(i q (Subscript[R, j,A])))}, {(Subscript[u, B])(e^(i q (Subscript[R, j,B])))}})),
=-(ω^2)(e^(-i ω t))(∑_j)(({{(Subscript[m, A])(Subscript[δ, i,j]), 0}, {0, (Subscript[m, B])(Subscript[δ, i,j])}}))(({{(Subscript[u, A])(e^(i q (2j)))}, {(Subscript[u, B])(e^(i q (2j+1)))}}))
=-(ω^2)(e^(-i ω t))(∑_j)(({{(Subscript[m, A])(Subscript[δ, i,j])(e^(i q (2j))), 0}, {0, (Subscript[m, B])(Subscript[δ, i,j])(e^(i q (2j+1)))}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}}))
=-(ω^2)(e^(-i ω t))(({{(Subscript[m, A])(e^(i q (2i))), 0}, {0, (Subscript[m, B])(e^(i q (2i+1)))}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}}))
=-(ω^2)(e^(-i ω t+2i q i))(({{(Subscript[m, A]), 0}, {0, (Subscript[m, B])(e^(i q))}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}})),

rhs=-(∑_j)(Subscript[K, i,j])(Subscript[x, j])=-(∑_j)(Subscript[K, i,j])(e^(-i ω t))(({{(Subscript[u, A])(e^(i q (Subscript[R, j,A])))}, {(Subscript[u, B])(e^(i q (Subscript[R, j,B])))}}))
=-κ (e^(-i ω t))(∑_j)(({{2(Subscript[δ, i,j]), -(Subscript[δ, i,j])-(Subscript[δ, i,j+1])}, {-(Subscript[δ, i,j])-(Subscript[δ, i,j-1]), 2(Subscript[δ, i,j])}}))(({{(Subscript[u, A])(e^(i q (2j)))}, {(Subscript[u, B])(e^(i q (2j+1)))}}))
=-κ (e^(-i ω t))(∑_j)(({{2(Subscript[δ, i,j])(e^(i q (2j))), (-(Subscript[δ, i,j])-(Subscript[δ, i,j+1]))(e^(i q (2j+1)))}, {(-(Subscript[δ, i,j])-(Subscript[δ, i,j-1]))(e^(i q (2j))), 2(Subscript[δ, i,j])(e^(i q (2j+1)))}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}}))
=-κ (e^(-i ω t))(({{2(e^(i q (2i))), -(e^(i q (2i+1)))-(e^(i q (2i-1)))}, {-(e^(i q (2i)))-(e^(i q (2i+2))), 2(e^(i q (2i+1)))}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}}))
=-κ (e^(-i ω t+2i q i))(({{2, -(e^(i q ))-(e^(-i q ))}, {-1-(e^(2i q )), 2(e^(i q ))}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}})),

(ω^2)(({{(Subscript[m, A]), 0}, {0, (Subscript[m, B])(e^(i q))}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}}))=κ (({{2, -(e^(i q ))-(e^(-i q ))}, {-1-(e^(2i q )), 2(e^(i q ))}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}})).

(ω^2)(({{(Subscript[m, A]), 0}, {0, (Subscript[m, B])}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}}))=κ (({{2, -(e^(i q ))-(e^(-i q ))}, {-(e^(-i q))-(e^(i q )), 2}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}}))
=2κ (({{1, -cos q}, {-cos q, 1}}))(({{(Subscript[u, A])}, {(Subscript[u, B])}})).

(Subscript[ω, q,±])=(Sqrt[(κ/((Subscript[m, A])(Subscript[m, B])))])((((Subscript[m, A])+(Subscript[m, B])±(Sqrt[(Subscript[m, A]^2)+(Subscript[m, B]^2)+2(Subscript[m, A])(Subscript[m, B])cos(2q)])))^(1/2))

Eigenvalues[{2κ{{1,-Cos[q]},{-Cos[q],1}},{{mA,0},{0,mB}}}]

Eigenvalues.

For κ=1, (Subscript[m, A])=1, (Subscript[m, B])=2:

Block[{κ=1,mA=1,mB=2},
Plot[{Sqrt[κ/(mA mB)(mA+mB-Sqrt[mA^2+mB^2+2mA mB Cos[2q]])],Sqrt[κ/(mA mB)(mA+mB+Sqrt[mA^2+mB^2+2mA mB Cos[2q]])]},{q,-π/2,π/2},PlotStyle->Black,Epilog->{{AbsoluteThickness[1],{Dashed,Line[{{0,0},{1/Sqrt[2κ/(mA+mB)],1}}],Text["\!\(v\_0 q\)",0.5{1/Sqrt[2κ/(mA+mB)],1},{0,0.8},{1,1}]},{Red,Line[{{0,0},{0.05,1.8}}],Text["\!\(c q\)",{0,1},{-1.5,0}]}}},FrameLabel->{"\!\( q\)","\!\(ω\_q\)"},FrameTicks->{{Automatic,Automatic},{{π#1,If[Mod[2#,1]==0,π Rationalize[#1]],#3}&@@@Charting`ScaledTicks[{Identity,Identity}][-1/2,1/2],Automatic}},ImageSize->180]]

Acoustic phonon: the lower branch (Subscript[ω, q,-]) phonon modes (sound waves). 
Definition: the branch that exhibits a linear dispersion as q→0 ((Subscript[ω, q->0]) is vanishing).

Derive Eq. () from Eq. ().

cos(2q)=1-((((2q))^2)/2)+𝒪((q^4)).

(Subscript[ω, q,-])=(Sqrt[(κ/((Subscript[m, A])(Subscript[m, B])))])((((Subscript[m, A])+(Subscript[m, B])-(Sqrt[(Subscript[m, A]^2)+(Subscript[m, B]^2)+2(Subscript[m, A])(Subscript[m, B])(1-((((2q))^2)/2))])))^(1/2))
=(Sqrt[(κ/((Subscript[m, A])(Subscript[m, B])))])((((Subscript[m, A])+(Subscript[m, B])-(Sqrt[(Subscript[m, A]^2)+(Subscript[m, B]^2)+2(Subscript[m, A])(Subscript[m, B])-4(Subscript[m, A])(Subscript[m, B])(q^2)])))^(1/2))
=(Sqrt[(κ/((Subscript[m, A])(Subscript[m, B])))])((((Subscript[m, A])+(Subscript[m, B])-((Subscript[m, A])+(Subscript[m, B]))(Sqrt[1-((4(Subscript[m, A])(Subscript[m, B]))/((((Subscript[m, A])+(Subscript[m, B])))^2))(q^2)])))^(1/2))
=(Sqrt[((κ ((Subscript[m, A])+(Subscript[m, B])))/((Subscript[m, A])(Subscript[m, B])))])(((1-(Sqrt[1-((4(Subscript[m, A])(Subscript[m, B]))/((((Subscript[m, A])+(Subscript[m, B])))^2))(q^2)])))^(1/2)).

(Subscript[ω, q,-])=(Sqrt[((κ ((Subscript[m, A])+(Subscript[m, B])))/((Subscript[m, A])(Subscript[m, B])))])(((1-(1-((2(Subscript[m, A])(Subscript[m, B]))/((((Subscript[m, A])+(Subscript[m, B])))^2))(q^2))))^(1/2))
=(Sqrt[((κ ((Subscript[m, A])+(Subscript[m, B])))/((Subscript[m, A])(Subscript[m, B])))])(((((2(Subscript[m, A])(Subscript[m, B]))/((((Subscript[m, A])+(Subscript[m, B])))^2))(q^2)))^(1/2))
=(Sqrt[((2κ)/((Subscript[m, A])+(Subscript[m, B])))])((((q^2)))^(1/2))
=(Sqrt[((2κ)/((Subscript[m, A])+(Subscript[m, B])))])|q|.

Optical phonon: the upper branch (Subscript[ω, q,+]) phonon modes (sound waves). 
Definition: the branch that maintains a finite frequency ((Subscript[ω, q->0]) is non-vanishing), such that the mode can possibly couples to light (photon). Because the speed of light c is much larger than the speed of sound (Subscript[v, 0]), such that the photon can transfer its energy-momentum to phonon only in the optical branch (hence the name “optical”).

There is a spectral gap between the two branches of sound waves.

(Subscript[x, i,A])=(Subscript[u, A])(e^(-i (Subscript[ω, q,±]) t+i q (Subscript[R, i,A]))), (Subscript[x, i,B])=(Subscript[u, B])(e^(-i (Subscript[ω, q,±]) t+i q (Subscript[R, i,B])))

Acoustic mode at q=0 (BZ center),

eigen value:(Subscript[ω, 0,-])=0,
eigen vector:(({{(Subscript[u, A])}, {(Subscript[u, B])}}))∝(({{1}, {1}})),

meaning that

(Subscript[x, i,A])=(Subscript[x, i,B])=u,

which corresponds to a mode that translate all atoms uniformly (the long-wave length limit of sound wave).

Acoustic mode at q=π/2 (BZ boundary), (assuming (Subscript[m, A])<(Subscript[m, B]))

eigen value:(Subscript[ω, π/2,-])=(Sqrt[((2κ)/(Subscript[m, B]))]),
eigen vector:(({{(Subscript[u, A])}, {(Subscript[u, B])}}))∝(({{0}, {1}})),

meaning that only B atoms oscillate

(Subscript[x, i,A])=0, (Subscript[x, i,B])=u (e^(-i (Subscript[ω, π/2,-]) t+i (π/2) (Subscript[R, i,B]))).

DynamicModule[{L=8,r=0.1,t=0,dt=0,x},
x[i_]:=If[EvenQ[i],0.3Cos[t+i π/2],0];
EventHandler[Dynamic[t+=dt;Graphics[{{LightGray,Line[{#,0}&/@{0.2,L+0.8}]},Translate[{{{EdgeForm[Black],Disk[{x[#],0},r(1+0.2(-1)^#)]}},{AbsoluteThickness[1],Line[{{0,-1.5r},{0,-2.5r}}],If[x[#]≠0,Line[{{{0,-2r},{x[#],-2r}},{{x[#]-0.5r Sign[x[#]],-1.5r},{x[#],-2r},{x[#]-0.5r Sign[x[#]],-2.5r}}}]]}},{#,0}]&/@Range[L]},ImageSize->300]],"MouseClicked":>If[dt>0,dt=0,dt=0.05π]]]

Optical mode at q=π/2 (BZ boundary), (assuming (Subscript[m, A])<(Subscript[m, B]))

eigen value:(Subscript[ω, π/2,+])=(Sqrt[((2κ)/(Subscript[m, A]))]),
eigen vector:(({{(Subscript[u, A])}, {(Subscript[u, B])}}))∝(({{1}, {0}})),

meaning that only A atoms oscillate

(Subscript[x, i,A])=u (e^(-i (Subscript[ω, π/2,+]) t+i (π/2) (Subscript[R, i,A]))),(Subscript[x, i,B])=0.

DynamicModule[{L=8,r=0.1,t=0,dt=0,x},
x[i_]:=If[OddQ[i],0.3Sin[t+i π/2],0];
EventHandler[Dynamic[t+=dt;Graphics[{{LightGray,Line[{#,0}&/@{0.2,L+0.8}]},Translate[{{{EdgeForm[Black],Disk[{x[#],0},r(1+0.2(-1)^#)]}},{AbsoluteThickness[1],Line[{{0,-1.5r},{0,-2.5r}}],If[x[#]≠0,Line[{{{0,-2r},{x[#],-2r}},{{x[#]-0.5r Sign[x[#]],-1.5r},{x[#],-2r},{x[#]-0.5r Sign[x[#]],-2.5r}}}]]}},{#,0}]&/@Range[L]},ImageSize->300]],"MouseClicked":>If[dt>0,dt=0,dt=0.05π Sqrt[2]]]]

At the BZ boundary, the oscillation mode transfers from B site (heavy atom, slower) to A site (light atom, faster), leading to a jump in the oscillation frequency => explains the spectral gap.

Optical mode at q=0 (BZ center),

eigen value:(Subscript[ω, 0,+])=(Sqrt[2κ((1/(Subscript[m, A]))+(1/(Subscript[m, B])))]),
eigen vector:(({{(Subscript[u, A])}, {(Subscript[u, B])}}))∝(({{(Subscript[m, B])}, {-(Subscript[m, A])}})),

meaning that

(Subscript[x, i,A])=(Subscript[m, B])u (e^(-i (Subscript[ω, 0,+]) t)),(Subscript[x, i,B])=-(Subscript[m, A])u (e^(-i (Subscript[ω, 0,+]) t)).

DynamicModule[{L=8,r=0.1,t=0,dt=0,x},
x[i_]:=If[OddQ[i],0.3Cos[t],-0.15Cos[t]];
EventHandler[Dynamic[t+=dt;Graphics[{{LightGray,Line[{#,0}&/@{0.2,L+0.8}]},Translate[{{{EdgeForm[Black],Disk[{x[#],0},r(1+0.2(-1)^#)]}},{AbsoluteThickness[1],Line[{{0,-1.5r},{0,-2.5r}}],If[x[#]≠0,Line[{{{0,-2r},{x[#],-2r}},{{x[#]-0.5r Sign[x[#]],-1.5r},{x[#],-2r},{x[#]-0.5r Sign[x[#]],-2.5r}}}]]}},{#,0}]&/@Range[L]},ImageSize->300]],"MouseClicked":>If[dt>0,dt=0,dt=0.05π Sqrt[3]]]]

Eigendecomposition of a matrix. Wikipedia.

Longitudinal and Transverse Modes

(Subscript[R, n])=(Subscript[n, 1])(Subscript[a, 1])+(Subscript[n, 2])(Subscript[a, 2])+(Subscript[n, 3])(Subscript[a, 3]),
({{{, ({{(Subscript[a, 1])=(1,0,0)}, {(Subscript[a, 2])=(0,1,0)}, {(Subscript[a, 3])=(0,0,1)}})}}).

H=(∑_i)((Subscript[p, i]^2)/(2m))+(1/2)(∑_(⟨i,j⟩))((Subscript[κ, ⊥]) ((((Subscript[x, i])-(Subscript[x, j])))^2)+((Subscript[κ, ∥])-(Subscript[κ, ⊥]))((((Subscript[(OverHat[R]), i,j])·((Subscript[x, i])-(Subscript[x, j]))))^2) ).

A lattice site is labeled by the index i (or j), which is equivalent to the integer vector n∈(ℤ^3). ⟨i,j⟩ denotes that (i,j) are across the nearest neighboring bond on the cubic lattice.

(Subscript[κ, ∥]) - parallel elastic constants of the bond (Young’s modulus),
(Subscript[κ, ⊥]) - perpendicular elastic constants of the bond (shear modulus),
such that suppose the bond (Subscript[(OverHat[R]), i,j]) is along x-direction, the bond potential energy is given by

(1/2)((Subscript[x, i])-(Subscript[x, j]))·(({{(Subscript[κ, ∥]),  ,  }, { , (Subscript[κ, ⊥]),  }, { ,  , (Subscript[κ, ⊥])}}))·((Subscript[x, i])-(Subscript[x, j])),

similar for y- and z-directions (by cyclic permuting the diagonal).

(Subscript[K, i,j])=(Subscript[K, 0]) (Subscript[δ, (Subscript[R, i]),(Subscript[R, j])])-(∑_(l=1,2,3))(Subscript[K, l]) ((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])])),
(Subscript[K, 0])=2((Subscript[κ, ∥])+2(Subscript[κ, ⊥]))𝟙 ,
(Subscript[K, l])=(Subscript[κ, ⊥])𝟙+((Subscript[κ, ∥])-(Subscript[κ, ⊥]))(Subscript[a, l])⊗(Subscript[a, l]),

(Subscript[K, 0])=2((Subscript[κ, ∥])+2(Subscript[κ, ⊥]))(({{1,  ,  }, { , 1,  }, { ,  , 1}})),
(Subscript[K, 1])=(({{(Subscript[κ, ∥]),  ,  }, { , (Subscript[κ, ⊥]),  }, { ,  , (Subscript[κ, ⊥])}})),
(Subscript[K, 2])=(({{(Subscript[κ, ⊥]),  ,  }, { , (Subscript[κ, ∥]),  }, { ,  , (Subscript[κ, ⊥])}})),
(Subscript[K, 3])=(({{(Subscript[κ, ⊥]),  ,  }, { , (Subscript[κ, ⊥]),  }, { ,  , (Subscript[κ, ∥])}})).

Derive Eq. () from Eq. ().

(∑_(⟨i,j⟩))((Subscript[κ, ⊥])(Subscript[x, i,j]^2)+((Subscript[κ, ∥])-(Subscript[κ, ⊥]))((((Subscript[(OverHat[R]), i,j])·(Subscript[x, i,j])))^2) )
=(1/2)(∑_(i,j))(∑_(l=1,2,3))((Subscript[κ, ⊥]) (Subscript[x, i,j]^2)+((Subscript[κ, ∥])-(Subscript[κ, ⊥]))((((Subscript[(OverHat[R]), i,j])·(Subscript[x, i,j])))^2) )((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])]))
=(1/2)(∑_(i,j))(∑_(l=1,2,3))((Subscript[κ, ⊥]) (Subscript[x, i,j]^2)+((Subscript[κ, ∥])-(Subscript[κ, ⊥]))((((Subscript[a, l])·(Subscript[x, i,j])))^2) )((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])]))
=(1/2)(∑_(i,j))(∑_(l=1,2,3))(Subscript[x, i,j])·((Subscript[κ, ⊥]) 𝟙+((Subscript[κ, ∥])-(Subscript[κ, ⊥]))(Subscript[a, l])⊗(Subscript[a, l]) )·(Subscript[x, i,j]) ((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])]))

=(1/2)(∑_(i,j))(∑_(l=1,2,3))((Subscript[x, i])-(Subscript[x, j]))·(Subscript[K, l])·((Subscript[x, i])-(Subscript[x, j])) ((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])]))
=(1/2)(∑_(i,j))(∑_(l=1,2,3))((Subscript[x, i])·(Subscript[K, l])·(Subscript[x, i])-(Subscript[x, i])·(Subscript[K, l])·(Subscript[x, j])-(Subscript[x, j])·(Subscript[K, l])·(Subscript[x, i])+(Subscript[x, j])·(Subscript[K, l])·(Subscript[x, j])) ((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])]))
=(∑_(i,j))(∑_(l=1,2,3))((Subscript[x, i])·(Subscript[K, l])·(Subscript[x, i])-(Subscript[x, i])·(Subscript[K, l])·(Subscript[x, j])) ((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])]))
=(∑_i)(∑_(l=1,2,3))2(Subscript[x, i])·(Subscript[K, l])·(Subscript[x, i])-(∑_(i,j))(∑_(l=1,2,3))(Subscript[x, i])·(Subscript[K, l])·(Subscript[x, j]) ((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])]))

=(∑_i)(Subscript[x, i])·(Subscript[K, 0])·(Subscript[x, i])-(∑_(i,j))(∑_(l=1,2,3))(Subscript[x, i])·(Subscript[K, l])·(Subscript[x, j]) ((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])]))
=(∑_(i,j))(Subscript[x, i])·(Subscript[K, 0])·(Subscript[x, j])(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])])-(∑_(i,j))(∑_(l=1,2,3))(Subscript[x, i])·(Subscript[K, l])·(Subscript[x, j]) ((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])]))
=(∑_(i,j))(Subscript[x, i])·((Subscript[K, 0])(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])])-(∑_(l=1,2,3))(Subscript[K, l])((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])])))·(Subscript[x, j]).

(Subscript[K, i,j])=(Subscript[K, 0])(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])])-(∑_(l=1,2,3))(Subscript[K, l])((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])])),

(Subscript[K, 0])=2(∑_(l=1,2,3))(Subscript[K, l])
=2(∑_(l=1,2,3))((Subscript[κ, ⊥]) 𝟙+((Subscript[κ, ∥])-(Subscript[κ, ⊥]))(Subscript[a, l])⊗(Subscript[a, l]))
=2(3(Subscript[κ, ⊥]) 𝟙+((Subscript[κ, ∥])-(Subscript[κ, ⊥]))(∑_(l=1,2,3))(Subscript[a, l])⊗(Subscript[a, l]))
=2(3(Subscript[κ, ⊥]) 𝟙+((Subscript[κ, ∥])-(Subscript[κ, ⊥]))𝟙)
=2((Subscript[κ, ∥])+2(Subscript[κ, ⊥]))𝟙.

Equation of motion:

m (Subscript[(Overscript[x, ..]), i])=-(∑_j)(Subscript[K, i,j])·(Subscript[x, j]).

Take the trial solution

(Subscript[x, i])=u (e^(-i ω t+i q·(Subscript[R, i]))),

substitute into Eq. (),

m (ω^2)u =(∑_j)(Subscript[K, i,j])·u (e^(-i q·((Subscript[R, i])-(Subscript[R, j]))))
=(∑_j)((Subscript[K, 0]) (Subscript[δ, (Subscript[R, i]),(Subscript[R, j])])-(∑_(l=1,2,3))(Subscript[K, l]) ((Subscript[δ, (Subscript[R, i]),(Subscript[R, j])+(Subscript[a, l])])+(Subscript[δ, (Subscript[R, i]),(Subscript[R, j])-(Subscript[a, l])])))·u (e^(-i q·((Subscript[R, i])-(Subscript[R, j]))))
=((Subscript[K, 0])-(∑_(l=1,2,3))(Subscript[K, l]) ((e^(-i q·(Subscript[a, l])))+(e^(i q·(Subscript[a, l])))))·u 
=((Subscript[K, 0])-(∑_(l=1,2,3))2(Subscript[K, l]) cos(q·(Subscript[a, l])))·u 
=K(q)·u,

given (Subscript[K, 0]) and (Subscript[K, l]) matrices in Eq. (),

K(q)=(({{(Subscript[K, x])(q),  ,  }, { , (Subscript[K, y])(q),  }, { ,  , (Subscript[K, z])(q)}})),
(Subscript[K, x])(q)=(Subscript[κ, ∥]) (((2sin ((Subscript[q, x])/2)))^2)+(Subscript[κ, ⊥]) (((2sin ((Subscript[q, y])/2)))^2)+(Subscript[κ, ⊥]) (((2sin ((Subscript[q, z])/2)))^2),
(Subscript[K, y])(q)=(Subscript[κ, ⊥]) (((2sin ((Subscript[q, x])/2)))^2)+(Subscript[κ, ∥]) (((2sin ((Subscript[q, y])/2)))^2)+(Subscript[κ, ⊥]) (((2sin ((Subscript[q, z])/2)))^2),
(Subscript[K, z])(q)=(Subscript[κ, ⊥]) (((2sin ((Subscript[q, x])/2)))^2)+(Subscript[κ, ⊥]) (((2sin ((Subscript[q, y])/2)))^2)+(Subscript[κ, ∥]) (((2sin ((Subscript[q, z])/2)))^2),

α labels the polarization of the phonon modes (sound waves) in the 3D solid.

Dispersion relation

Block[{κ1=1,κ2=0.5,ωs},
ωs[{qx_,qy_,qz_}]:=Sqrt[{{κ1,κ2,κ2},{κ2,κ1,κ2},{κ2,κ2,κ1}}.(2{Sin[qx/2],Sin[qy/2],Sin[qz/2]})^2];
BZPlot[ωs,{"Γ","X","M","Γ","R","X"},BZRules->{"Γ"->{0,0,0},"X"->{π,0,0},"M"->{π,π,0},"R"->{π,π,π}},Epilog->{{#3,Text["\!\(ω\_\(q,"<>#1<>"\)\)",#2]}&@@@{{"x",{0.5π,2.3},Blue},{"y",{2.8π,2.5},Red},{"z",{2.2π,1.5},Green}}},FrameLabel->{None,"\!\(ω\_\(q\)\)"},AspectRatio->0.4,ImageSize->240]]

plotted along high symmetry lines in the cubic lattice Brillouin zone

Graphics3D[{{FaceForm[Opacity[0.2]],EdgeForm[AbsoluteThickness[1]],Cube[{0,0,0},2]},{Black,Rotate[Arrow[Tube[{{0,0,0},{2,0,0}}]],2π #/3,{1,1,1}]&/@{0,1,2},
Text["\!\(b\_"<>#1<>"\)",##2]&@@@{{"1",{2,0,0},{0.1,0.3}},{"2",{0,2,0},{1,0.8}},{"3",{0,0,2},{1.2,1}}}},
GraphicsComplex[{{0,0,0},{1,0,0},{1,1,0},{1,1,1}},{{AbsolutePointSize[6],Point[{1,2,3,4}]},{AbsoluteThickness[2],Line[{1,2,3,1,4,2}],Line[{4,3}]},{Text["Γ",1,{1,-0.5}],Text["X",2,{-0.5,1.2}],Text["M",3,{-1.5,0.2}],Text["R",4,{0.5,-0.8}]}}]},Boxed->False,ViewPoint->{5,2,3},Lighting->{{"Directional",White,{5,0,0}},{"Directional",White,{0,5,5}}},ImageSize->160]

Longitudinal phonon: polarization parallel to momentum (u∥q)

Transverse phonon: polarization perpendicular to momentum (u⊥q)

Phonon modes along Γ-X, i.e. q=q(1,0,0),

Longitudinal mode (one branch)

(Subscript[ω, q,x])=2(Sqrt[(((Subscript[κ, ∥]) )/m)])|sin (q/2)|(Overscript[→, q->0])(Subscript[v, L])|q|.

Longitudinal speed of sound

(Subscript[v, L])=(Sqrt[(((Subscript[κ, ∥]) )/m)]).

Transverse mode (two branches)

(Subscript[ω, q,y])=(Subscript[ω, q,z])=2(Sqrt[(((Subscript[κ, ⊥]) )/m)])|sin (q/2)|(Overscript[→, q->0])(Subscript[v, T])|q|.

Transverse speed of sound

(Subscript[v, T])=(Sqrt[(((Subscript[κ, ⊥]) )/m)]).

Phonon modes along Γ-M, i.e. q=q(1,1,0)/(Sqrt[2]),

Longitudinal mode (one branch)

(Subscript[ω, q,x+y])=2(Sqrt[(((Subscript[κ, ∥])+(Subscript[κ, ⊥]) )/(2m))])|sin (q/2)|(Overscript[→, q->0])(Subscript[v, L])|q|.

Longitudinal speed of sound

(Subscript[v, L])=(Sqrt[(((Subscript[κ, ∥])+(Subscript[κ, ⊥]) )/(2m))]).

Transverse mode (two branches)

(Subscript[ω, q,x-y])=2(Sqrt[(((Subscript[κ, ∥])+(Subscript[κ, ⊥]) )/(2m))])|sin (q/2)|(Overscript[→, q->0])(Subscript[v, T1])|q|,
(Subscript[ω, q,z])=2(Sqrt[(((Subscript[κ, ⊥]) )/m)])|sin (q/2)|(Overscript[→, q->0])(Subscript[v, T2])|q|.

Transverse speed of sound

(Subscript[v, T1])=(Sqrt[(((Subscript[κ, ∥])+(Subscript[κ, ⊥]) )/(2m))]),
(Subscript[v, T2])=(Sqrt[(((Subscript[κ, ⊥]) )/m)]).

The degeneracy between (Subscript[v, L]) and (Subscript[v, T1]) is accidental. Including further neighboring bonds in the model will lift this degeneracy.

Specific Heat of Solid

Boltzmann Model

(Subscript[c, V])=((∂⟨ϵ⟩)/(∂T)).

Specific heat often varies with temperature. Its temperature dependence reflects how energy is stored as excitations in the material.

Energy of a harmonic oscillator

ϵ(x,p)=(1/(2m))(p^2)+(κ/2)(x^2).

The oscillation frequency is ω=(Sqrt[κ/m]).

Probability for the oscillator to be at the state (x,p) is

p(x,p)=(1/Z)(e^(-β ϵ(x,p))),

where Z is the partition function (classical, canonical ensemble)

Z=∫(d^D)p(d^D)x (e^(-β ϵ(x,p)))
=(((((2π m)/β)))^(D/2))(((((2π)/(β κ))))^(D/2))
=(((((2π)/(β ω))))^D).

β=1/((Subscript[k, B])T): (Subscript[k, B]) - Boltzmann constant, T - temperature.

D - dimension of space.

Average energy

⟨ϵ⟩=∫(d^D)p(d^D)xϵ(x,p)p(x,p)
=-(∂_β) ln Z
=D (Subscript[k, B])T.

Specific heat

For 3D materials (D=3),

(Subscript[c, V])=3(Subscript[k, B]),

which is known as the Dulong-Petit law. Specific heats of some solids at room temperature and pressure:

({{Material, (Subscript[c, V])/(Subscript[k, B]),  }, {Gold (Au), 3.05,  }, {Antimony (Sb), 3.03,  }, {Zinc (Zn), 3.02,  }, {Iron (Fe), 3.02,  }, {Silver (Ag), 2.99,  }, {Lithium (Li), 2.98,  }, {Copper (Cu), 2.94,  }, {Aluminum (Al), 2.91,  }, {Diamond (C), 0.735, (why?)}})

Einstein Model

Energy is quantized in quantum mechanics. For D-dimensional harmonic oscillator, its eigen energies (energy levels) are labeled by n∈(ℕ^D),

(Subscript[ϵ, n])=ℏ ω (∑_(i=1)^D)((Subscript[n, i])+(1/2)).

ω - oscillator frequency,

ℏ - Planck constant.

Probability to find the oscillator in the (Subscript[ϵ, n]) energy level

(Subscript[p, n])=(1/Z)(e^(-β (Subscript[ϵ, n]))),

where Z is the partition function (quantum, canonical ensemble)

Z=(∑_(n∈(ℕ^D)))(e^(-β (Subscript[ϵ, n])))
=(((2sinh(β ℏ ω/2)))^(-D)).

Evaluate the partition function in Eq. ().

Z=(∑_(n∈(ℕ^D)))(e^(-β (Subscript[ϵ, n])))
=(∑_(n∈(ℕ^D)))exp(-β ℏ ω (∑_(i=1)^D)((Subscript[n, i])+(1/2)))
=(∑_(n∈(ℕ^D)))(∏_(i=1)^D)exp(-β ℏ ω ((Subscript[n, i])+(1/2)))
=(∏_(i=1)^D)(∑_((Subscript[n, i])∈ℕ))(e^(-β ℏ ω ((Subscript[n, i])+(1/2))))
=(∏_(i=1)^D)((e^(-(1/2)β ℏ ω ))(∑_((Subscript[n, i])=0)^∞)(e^(-β ℏ ω (Subscript[n, i]))))
=(∏_(i=1)^D)((e^(-(1/2)β ℏ ω ))(1/(1-(e^(-β ℏ ω)))))
=(∏_(i=1)^D)((1/(2sinh(β ℏ ω/2))))
=(((2sinh(β ℏ ω/2)))^(-D)).

Average energy

⟨ϵ⟩=(∑_(n∈(ℕ^D)))(Subscript[ϵ, n])(Subscript[p, n])
=-(∂_β) ln Z
=((D ℏ ω)/2)coth(((β ℏ ω)/2))
=D ℏ ω ((Subscript[n, B])(β ℏ ω)+(1/2)),

(Subscript[n, B]) is the Bose-Einstein distribution function

describing the average number of bosons occupying the mode of energy ℏ ω.

Specific heat

(Subscript[c, V])=((∂⟨ϵ⟩)/(∂T))
=D (Subscript[k, B]) ((((((β ℏ ω)/2))csch(((β ℏ ω)/2))))^2).

Define a temperature scale (the Einstein temperature)

(Subscript[T, Θ])=((ℏ ω)/(Subscript[k, B])),

Eq. () can be written as

Off[General::munfl];Plot[(1/(2x)Csch[1/(2x)])^2,{x,0,2},PlotStyle->Black,FrameLabel->{"\!\(T/T\_Θ\)","\!\(c\_V/D k\_B\)"},ImageSize->180]

In the high-temperature limit (T≫(Subscript[T, Θ])), the Dulong-Petit law (Subscript[c, V])=D (Subscript[k, B]) is recovered.

However, at low-temperature (T≪(Subscript[T, Θ])), the degrees of freedom “freeze out”. The energy level is discrete in quantum mechanics, when the temperature is too low compared to the excitation energy (the level spacing), the system will stuck in the ground state and can not absorb/release energy, hence no contribution to the specific heat.

Debye Model

Phonon can have different frequencies, as described by the dispersion relation (Subscript[ω, q,α]), which depends on:

q - phonon momentum,

α - phonon band index (including acoustic/optical branches, and different polarizations).

Following Eq. (), the specific heat should also average over all phonon modes

N - number of atoms, which is also the number of phonon modes (per polarization).

At low-temperature, focus on acoustic phonons, its dispersion can be modeled by

(Subscript[ω, q])=v|q|,

where v is the speed of sound.

Here the dispersion is assumed to be isotropic, ignoring the velocity difference between longitudinal and transverse phonons.

The frequency is bounded as q varies in the Brillouin zone. So there must be some maximal frequency (Subscript[ω, max])  (to be specified later).

Under the isotropic assumption, the momentum summation can be written as integration

(1/N)(∑_α)(∑_(q∈BZ))=(D/N)((((L/(2π))))^D)∫(d^D)q
=((D Ω)/(((2π))^D))∫(Subscript[A, D])(q^(D-1))dq
=∫((D (Subscript[A, D])Ω)/(((2π v))^D))(Subscript[ω, q]^(D-1))d(Subscript[ω, q])
=∫d(Subscript[ω, q]) g((Subscript[ω, q]))

D - space dimension (also the number of polarizations α=1,…,D)

L - linear size of the system, which determined the momentum discretization unit (2π/L).

Ω - volume of unit cell.

(Subscript[A, D])=2(π^(D/2))/Γ((D/2)) - area of a (D-1)-dimensional unit sphere (in the D-dimensional space).

Define the phonon density of state

where the Debye frequency (Subscript[ω, Θ]) is introduced as a frequency scale

(Subscript[ω, Θ]^D)=((D (((2π v))^D))/((Subscript[A, D])Ω)),

or more explicitly

(Subscript[ω, Θ])=({{{, ({{(π/Ω)v, D=1}, {(((((4π)/Ω)))^(1/2))v, D=2}, {(((((6(π^2))/Ω)))^(1/3))v, D=3}, {…, …}})}})

Table[𝔇->PowerExpand[(𝔇(2π v)^𝔇/(2π^(𝔇/2)/Gamma[𝔇/2] Ω))],{𝔇,3}]

Enumerate (Subscript[ω, Θ]).

The maximal frequency (Subscript[ω, max]) must be such that the density of state integrates to the number of polarization modes

D=(∫_0^(Subscript[ω, max]))dω g(ω)=D (((((Subscript[ω, max])/(Subscript[ω, Θ]))))^D),

so it turns out that the Debye frequency is also the maximal frequency (Subscript[ω, max])=(Subscript[ω, Θ]).

The specific heat Eq. () can then be evaluated as

(Subscript[c, V])=(Subscript[k, B]) (∫_0^(Subscript[ω, Θ]))dω g(ω)((((((β ℏ ω)/2))csch(((β ℏ ω)/2))))^2),

which can be organized as

(Subscript[T, Θ]) - Debye temperature, define as

(Subscript[T, Θ])=((ℏ (Subscript[ω, Θ]))/(Subscript[k, B])).

Block[{𝔇=3,Y,Cv},
Y[y_]:=𝔇 NIntegrate[x^(𝔇-1)(x/2 Csch[x/2])^2,{x,0,y}];
Cv[T_]:=T^𝔇 Y[1/T];
Plot[{Cv[T],(1/(2T)Csch[1/(2T)])^2},{T,0.001,1},PlotStyle->{Black,Directive[Black,Dashed]},FrameLabel->{"\!\(T/T\_Θ\)","\!\(c\_V/D k\_B\)"},Epilog->{Text["Debye",{0.22,0.8}],Text["Einstein",{0.57,0.4}]},ImageSize->180]]

Low-temperature limit T≪(Subscript[T, Θ]),

The specific heat Eq. () reads

(Subscript[c, V])=D (((Subscript[k, B])((T/(Subscript[T, Θ]))))^D)Y(∞)
=(((Subscript[k, B])((T/(Subscript[T, Θ]))))^D)(D^2) (D+1)!ζ(D+1)

where ζ is the Riemann zeta function, or more explicitly

(Subscript[c, V])=(Subscript[k, B])({{{, ({{((T/(Subscript[T, Θ])))((π^2)/3), D=1}, {((((T/(Subscript[T, Θ]))))^2)24ζ(3), D=2}, {((((T/(Subscript[T, Θ]))))^3)((12(π^4))/5), D=3}, {…, …}})}}).

Table[𝔇->𝔇^2 Gamma[2+𝔇] PolyLog[1+𝔇,1],{𝔇,3}]

Enumerate (Subscript[C, V]).

In conclusion, the low-temperature specific heat from the acoustic phonon contribution scales with temperature as

which is (Subscript[c, V])∼(T^3) for D=3.

Stability of Crystal Order

Boltzmann Model

Energy of a harmonic oscillator

ϵ(x,p)=(1/(2m))(p^2)+(κ/2)(x^2).

Probability for the oscillator to be at the state (x,p) is

p(x,p)=(1/Z)(e^(-β ϵ(x,p))),

where Z is the partition function

Z=∫(d^D)p(d^D)x (e^(-β ϵ(x,p)))
=(((((2π m)/β)))^(D/2))(((((2π)/(β κ))))^(D/2)).

β=(Subscript[k, B])T: (Subscript[k, B]) - Boltzmann constant, T - temperature.

D - dimension of space.

Variance of displacement

⟨(x^2)⟩=∫(d^D)p(d^D)x (x^2)p(x,p)
=-(2/β)(∂_κ)ln Z
=((D (Subscript[k, B])T)/κ).

The result is consistent with the equal partition theorem, that the average energy stored in each quadratic term of the Hamiltonian is (1/2)(Subscript[k, B ])T. Thus, for harmonic oscillator the average potential energy should be

⟨(1/2)κ (x^2)⟩=(D/2)(Subscript[k, B])T

Strictly speaking, the variance should be given by ⟨(x^2)⟩-((⟨x⟩)^2), but the first moment ⟨x⟩=0 vanishes due to the inversion symmetry ϵ(x,p)=ϵ(-x,p) of the harmonic oscillator.

Melting condition: If the displacement variance exceeds the squared atomic spacing (a^2) of the lattice, the solid will melt.

⟨(x^2)⟩=((D (Subscript[k, B])T)/κ)>(a^2)=>T>((κ (a^2))/(D (Subscript[k, B])))=(SubStar[T]).

The melting temperature (SubStar[T]) is consistent with the rough estimate in Eq. () (SubStar[T])=(Subscript[E, bind])/(N (Subscript[k, B])), as (Subscript[E, bind])∼N κ (a^2).

However, this estimation did not takes into account the quantum effect and the collective nature of phonon excitations, and is wrong in low dimensions.

Einstein Model

Eigen energies of a D-dimensional harmonic oscillator (labeled by n∈(ℕ^D))

(Subscript[ϵ, n])=ℏ ω (∑_(i=1)^D)((Subscript[n, i])+(1/2)).

ω=(Sqrt[κ/m]) - oscillator frequency.

ℏ - Planck constant.

Probability to find the oscillator in the (Subscript[ϵ, n]) energy level

(Subscript[p, n])=(1/Z)(e^(-β (Subscript[ϵ, n]))),

where Z is the partition function (quantum, canonical ensemble)

Z=(∑_(n∈(ℕ^D)))(e^(-β (Subscript[ϵ, n])))
=(((2sinh(β ℏ ω/2)))^(-D)).

Variance of displacement on a particular eigenstate (Ket[n])

(Bra[n])(x^2)(Ket[n])=2(∂_κ)(Subscript[ϵ, n]).

Variance of displacement in the thermal ensemble

⟨(x^2)⟩=(∑_(n∈(ℕ^D)))(Bra[n])(x^2)(Ket[n])(Subscript[p, n])
=(∑_(n∈(ℕ^D)))2((∂_κ)(Subscript[ϵ, n]))(Subscript[p, n])
=-(2/β)(∂_κ)ln Z
=((D ℏ)/(2m ω))coth(((β ℏ ω)/2))
=((D ℏ ω)/κ)((Subscript[n, B])(β ℏ ω)+(1/2)).

Introduce the Einstein temperature following Eq. ()

(Subscript[T, Θ])=((ℏ ω)/(Subscript[k, B])),

the variance Eq. () can be written as

Plot[{1/(Exp[1/x]-1)+1/2,x},{x,0,2},PlotStyle->{Black,Directive[Black,Dashed]},FrameLabel->{"\!\(T/T\_Θ\)","\!\(⟨x\^2⟩/x\_Θ\%2\)"},Epilog->{Text["Eistein",{0.5,1.2}],Text["Boltzmann",{1.4,1.4},{0,1},{1,0.618}]},ImageSize->180]

At high temperature (T≫(Subscript[T, Θ])), the linear temperature behavior is restored,

⟨(x^2)⟩=(Subscript[x, Θ]^2)((T/(Subscript[T, Θ])))=((D (Subscript[k, B])T)/κ),

matching the Boltzmann model Eq. ().

At low temperature (T≪(Subscript[T, Θ])), the fluctuation saturates to a finite level

⟨(x^2)⟩=((Subscript[x, Θ]^2)/2)=((D (Subscript[k, B])(Subscript[T, Θ]))/(2κ))=((D ℏ)/(2m ω)).

This non-vanishing fluctuation at zero temperature is called the zero-point fluctuation, or the quantum fluctuation, originated from the uncertainty relation between position and momentum in quantum mechanics.

Debye Model

Define the average displacement variance of an atom in the solid

⟨(x^2)⟩=(1/N)(∑_i)⟨(Subscript[x, i]^2)⟩.

Fourier transform to the momentum space

(Subscript[x, i])=(1/(N^(1/2)))(∑_(q∈BZ))(Subscript[x, q])(e^(i q·(Subscript[R, i])))

The fluctuation Eq. () can be written as

⟨(x^2)⟩=(1/N)(∑_(q∈BZ))⟨(Subscript[x, q]^†)(Subscript[x, q])⟩.

Each momentum q labels an independent acoustic phonon mode (for each polarization), for which (Subscript[x, q]) fluctuates as an independent harmonic oscillator.

Use the quantum result in Eq. (), and take the Debye model for density of state in Eq. (),

⟨(x^2)⟩=(1/N)(∑_(q∈BZ))((D ℏ)/(2m (Subscript[ω, q])))coth(((β ℏ (Subscript[ω, q]))/2))
=(∫_0^(Subscript[ω, Θ]))dω g(ω)(ℏ/(2m ω))coth(((β ℏ ω)/2)),

which can be organized as

For D>2, at high temperature (T≫(Subscript[T, Θ])), ⟨(x^2)⟩ approaches to the linear T behavior

⟨(x^2)⟩=(Subscript[x, Θ]^2)((T/(Subscript[T, Θ])))=(D/((D-2)m (Subscript[ω, Θ]^2)))D (Subscript[k, B])T.

Block[{𝔇=3,Z},
Z[y_]:=NIntegrate[x^(𝔇-2)/2Coth[x/2],{x,0.y,y}];
Plot[{(𝔇-2)x^(𝔇-1) Z[1/x],x},{x,0.001,2},PlotStyle->{Black,Directive[Black,Dashed]},FrameLabel->{"\!\(T/T\_Θ\)","\!\(⟨x\^2⟩/x\_Θ\%2\)"},Epilog->{Text["Debye",{0.5,1.}],Text["Boltzmann",{1.4,1.4},{0,1},{1,0.618}]},PlotLabel->StringTemplate["\!\(D=`1`\)"][𝔇],ImageSize->180]]

In the zero temperature limit (T→0), the fluctuation saturates to

⟨(x^2)⟩=((D-2)/(2(D-1)))(Subscript[x, Θ]^2)=(((D^2) ℏ)/(2(D-1)m (Subscript[ω, Θ]))).

Derive Eq. () and Eq. ().

Z(((Subscript[T, Θ])/T)→0)=(D-2)(∫_0^((Subscript[T, Θ])/T))dx ((x^(D-2))/2)(2/x)
=(D-2)(∫_0^((Subscript[T, Θ])/T))dx (x^(D-3))
= (Subscript[((x^(D-2))), x=0]^(x=(Subscript[T, Θ])/T))
=(((((Subscript[T, Θ])/T)))^(D-2)),

⟨(x^2)⟩=(((Subscript[x, Θ]^2)((T/(Subscript[T, Θ]))))^(D-1))(((((Subscript[T, Θ])/T)))^(D-2))
=(Subscript[x, Θ]^2)((T/(Subscript[T, Θ]))).

Z(((Subscript[T, Θ])/T)→∞)=(D-2)(∫_0^((Subscript[T, Θ])/T))dx ((x^(D-2))/2)
=((D-2)/(2(D-1)))(Subscript[((x^(D-1))), x=0]^(x=(Subscript[T, Θ])/T))
=((D-2)/(2(D-1)))(((((Subscript[T, Θ])/T)))^(D-1)),

⟨(x^2)⟩=(((Subscript[x, Θ]^2)((T/(Subscript[T, Θ]))))^(D-1))((D-2)/(2(D-1)))(((((Subscript[T, Θ])/T)))^(D-1))
=((D-2)/(2(D-1)))(Subscript[x, Θ]^2).

For D=2, according to Eq. (), ⟨(x^2)⟩ diverges at any finite temperature => 2D crystals are inevitably melted by thermal fluctuations.
However, 2D crystals can still be stable at zero-temperature, as ⟨(x^2)⟩ remains finite in the T→0 limit for D=2 in Eq. ().

For D=1, according to Eq. (), ⟨(x^2)⟩ even diverges in the zero-temperature limit => 1D crystals are inevitably melted by quantum fluctuations.

Quantum-classical duality: a quantum system (T=0) at D dimension has properties similar to a classical system (T>0) at (D+1) dimension. This duality connects quantum mechanics with statistical mechanics in one higher dimension.