Filling a triangular lattice

(a) Reciprocal lattice vectors

First, identify our two crystal vectors.

a1={(3/2),-((Sqrt[3])/2)};
a2={(3/2),((Sqrt[3])/2)};

Go to 3D so that we can simply calculate the reciprocal lattice vectors.

a1temp=Append[a1,0];
a2temp=Append[a2,0];
a3temp={0,0,1};

vol=a3temp.Cross[a1temp,a2temp];
b1temp=2π Cross[a2temp,a3temp]/vol;
b2temp=2 π Cross[a3temp,a1temp]/vol;

Make them 2D again.   Our answer is in units of 1/a

b1=Take[b1temp,2]
b2=Take[b2temp,2]

(b)  Plotting k-space and Brillouin zones.

Make a list of points in K-space, and plot them.

kspace=Flatten[Table[m b1 + n b2,{m,-7,7},{n,-5,5}],1];
kPlot=ListPlot[kspace,AspectRatio-> 1,PlotRange-> {{-15,15},{-15,15}},PlotStyle-> PointSize[Large]]

ListPlot[kspace,AspectRatio-> 1,PlotStyle-> PointSize[Large],PlotRange-> {{-40,40},{-40,40}}]

Make the lines that define the BZ edge:  to do this use the point slope formula of a line.   We know a point on the BZ line (it's half the G vector) and the slope (it's perpendicular to the G vector).   This gives us enough to solve for the slope-intercept form of the BZ line.   This gives us "y = m x + b " where m is the slope and b is the intercept.  We have to be careful with vertical lines which have infinite slope.   We also have to drop the line generated by the RL vector G=0.

lines=Chop[Flatten[Table[ pt= (m b1 + n b2)/2 ;  
slope = - pt[[1]]/(pt[[2]]+.00001); 
b=pt[[2]] - slope * pt[[1]];
slope * x + b,{m,-7,7},{n,-7,7}],1]];
zeropos=Position[lines,0][[1]];
bzEdges=Drop[lines,zeropos];
bzPlot=Plot[bzEdges,{x,-10,10},AspectRatio-> 1,PlotRange-> {{-10,10},{-10,10}},PlotPoints-> Automatic]

Show[kPlot,bzPlot]

Ok, here's a cool color map.   We make it by creating a function that jumps by one every time we cross a RL line.  For any x and y we calculate where a line drawn from the origin to (x,y) crosses the BZ line. We then calculate if the distance to the point of intersection is less than the distance to (x,y).  If it is, then (x,y) is beyond this BZ edge.  We then count how many BZ edges we cross.

fxns=Chop[Flatten[Table[ pt= (m b1 + n b2)/2 ;  
slope = - pt[[1]]/(pt[[2]]+.00001); 
b=pt[[2]] - slope * pt[[1]];
(b/(y-slope x)),{m,-5,5},{n,-5,5}],1]];
zeropos=Position[fxns,0][[1]];
stepArgs=Drop[fxns,zeropos];
temp8=Sum[UnitStep[((stepArgs[[i]])^2)-1],{i,1,Length[stepArgs]}];

rb=DensityPlot[temp8,{x,-7,7},{y,-6.5,6.5},ColorFunction-> "Rainbow",PlotPoints-> 64]

Plot this with the BZ lines above.   We can see that we have color coded the BZ's.   It's almost perfect.

Show[rb,bzPlot]

Plot3D[temp8,{x,-6.5,6.5},{y,-6.25,6.25},ColorFunction-> "Rainbow",PlotPoints-> 64]

(c)  Fermi surface intersections

We can determine the fermi wavevector by counting states.  Including the factor of two for spin:

	N = 2 (( π (kf^2))/(((2π/L))^2))
	
	n =  ((kf^( 2))/(2 π))
	
The volume of a cell in real space is 3 (a^2)/2.  If we have q electrons per unit cell then:

	kf = (Sqrt[2 π  (2/(3 (a^2))) s ])=(Sqrt[((4 π  s)/3)]) (1/a) 
	
Our kspace plots are in units of 1/a.   Let's draw the circles for s=1, 2 and 3:

fs=Show[Graphics[Circle[{0,0},(Sqrt[((4 π)/3)])]],Graphics[Circle[{0,0},(Sqrt[((8 π)/3)])]],Graphics[Circle[{0,0},(Sqrt[((12 π)/3)])]]]

Show[rb,fs]

The monovalent case is entirely within the first BZ,  the divalent is in the first and second, as is the trivalent case.