Code for article:
Fermi Surfaces Coloured with the Group Speed to Reveal Critical Points for Singularities of the Density of States

by
E.G. Emberly and P.D. Loly
Department of Physics
University of Manitoba
Winnipeg, Manitoba
Canada  R3T 2N2
loly@cc.umanitoba.ca

Published in
Mathematica in Education and Research
Vol.4 No.1

Copyright 1995 TELOS/Springer-Verlag

The first three commands of this code simply define the necessary equations for generating the fermi surface, with the last command defining a list of the desired energies to be plotted.  In order for the surface to be colored with a separate function (i.e. the group speed), a Plot3D command must be used.  Because of this, the equation for the fermi surface was solved for z, which yielded an inverse cosine.  Mathematica only provides one solution for ArcCos, so only the top of the surface can be plotted.

The following defines the surface function, which depends on kx, ky and also the energy, n.

f[n_, kx_, ky_]:=
	ArcCos[3 - n - Cos[kx] - Cos[ky]]

This command defines the colour function, with which the surface is to be coloured.

colour[n_,kx_,ky_]:=
	(Sin[kx]^2 + Sin[ky]^2 + Sin[ f[n, kx, ky] ]^2)^(1/2)

Here is the list of the energies to be plotted.

energy={1.5,2,2.5,3.5,4,4.5};

The next section of code generates the graphics array of the six previously defined energies.  A Table command is used to generate a list of six plots, one for each energy in the list.  In regards to the Hue command, it only takes an arguement which ranges between 0 and 1 (with 0 being red, then going to blue at around .75, and then back to red at 1).  It was desired to have the zeros of the group speed being represented by blue and the maximum value of the group speed being represented by red.  The function within the Hue command achieves this.
The options that follow are fairly standard, except for perhaps ClipFill.  This option causes Mathematica not to render any part of the surface which goes outside the bounding box, or any imaginary result that might be produced by the equation for the surface.  By using Partition the graphs are organized into two columns.  These are then displayed in a GraphicsArray with the Show command.  It was desired to achieve high resolution output for each of the surfaces.  This however takes a lot of computation time, and for everyday use, a lower PlotPoint specification may be more practical.

graphs1=Table[
	Plot3D[
		{f[energy[[i]], kx, ky],
		Hue[.75(1 - colour[energy[[i]], kx, ky]/Sqrt[3])]},
		{kx, -Pi, Pi},{ky, -Pi, Pi},
		PlotRange -> {{-Pi,Pi},{-Pi,Pi},{0,Pi}},
		Mesh -> False,
		PlotPoints -> 250,
		ClipFill -> None,
		DisplayFunction -> Identity],
	{i, 1, 6}];

$DefaultFont={"Helvetica",7};

Show[GraphicsArray[ Partition[graphs1,2] ],
	DisplayFunction -> $DisplayFunction]

fermi3 code

f[energy_,kx_,ky_]:=ArcCos[3-energy-Cos[kx]-Cos[ky]]

colour[energy_,kx_,ky_]:=
	Sqrt[Sin[kx]^2+Sin[ky]^2+Sin[f[energy,kx,ky]]^2]

Plot3D[{f[3,kx,ky],
	Hue[.75(1-colour[3,kx,ky]/Sqrt[3])]},
	{kx,-Pi,Pi},{ky,-Pi,Pi},
	ClipFill -> None,
	Mesh -> False,
	PlotPoints -> 180];