# Plotting 2-D Electron Density Carbon Atoms

I need to create a 2-D electron energy density plot in Mathematica to compare with my STM experimental results in my lab class. This would be done by plotting the sum of the symmetric and anti-symmetric wave functions,

$$\Psi_s(\textbf{r}) = \sum_{i=1}^{3}a_i\cos(\textbf{k}_i\cdot\textbf{r})$$ $$\Psi_a(\textbf{r}) = \sum_{i=1}^{3}b_i\sin(\textbf{k}_i\cdot\textbf{r})$$ and squaring it to get this, $$\Psi_s^2(\textbf{r})+\Psi_a^2(\textbf{r})=\Big[\sum_{i=1}^{3}a_i\cos(\textbf{k}_i\cdot\textbf{r})\Big]^2+\Big[\sum_{i=1}^{3}b_i\sin(\textbf{k}_i\cdot\textbf{r})\Big]^2$$ Where the k's are equal to, $$\textbf{k}_1=\frac{4\pi}{3d_b}\hat{y}$$ $$\textbf{k}_2=\frac{4\pi}{3d_b}(-\frac{1}{2}\sqrt{3}\hat{x}-\frac{1}{2}\hat{y})$$ $$\textbf{k}_3=\frac{4\pi}{3d_b}(\frac{1}{2}\sqrt{3}\hat{x}-\frac{1}{2}\hat{y})$$ Here is what I have so far

Needs["VectorAnalysis"]
r = {1, 1, 0};
k1 = 4*Pi/3 (2.46) {0, 1, 0};
k2 = 4*Pi/3 (2.46) {-(Sqrt[3]/2), -(1/2), 0};
k3 = 4*Pi/3 (2.46) {(Sqrt[3]/2), -(1/2), 0};

Subscript[\[Psi], s] =
Cos[DotProduct[k1, r]]*x + Cos[DotProduct[k2, r]]*x +
Cos[DotProduct[k3, r]]*x;

Subscript[\[Psi], a] =
Sin[DotProduct[k1, r]]*y + Sin[DotProduct[k2, r]]*y +
Sin[DotProduct[k3, r]]*y;
DensityPlot[(Subscript[\[Psi], s] x)^2 + (Subscript[\[Psi], a]
y)^2, {x, 0, 100}, {y, 0, 100}]


Resulting in this

However, it should be uniform distribution of "bright spots" to model the electron density of carbon atoms, something along the lines of this

Could someone help me out with this? It would be much appreciated.

• Hi Rosie, and welcome to Mathematica.SE! I hope you will become a regular contributor. Thank you for your complete and well-formatted post. I will post a complete answer shortly, but first, please note that you haven't actually specified functions of x and y, and that Subscript can get you into trouble. I strongly suggest that you have a look at the material here for some ideas. Also, what version are you using? The VectorAnalysis package is obsolete and deprecated. – Verbeia Sep 30 '15 at 3:53
• @Rosie, do not forget to tick (not only upvote) the answer if it solves your problem =)) – garej Sep 30 '15 at 20:24
• @garej Thanks I will do so – Rosie Sep 30 '15 at 21:06

If you define your constants this way, you can write the functions using Sum - and looking at your code, it seems that you wanted to define db=2.46 - if so, then you did not correctly put it in the denominator, it was in the numerator.

db = 2.46;
k[1] = 4*Pi/(3 db) {0, 1};
k[2] = 4*Pi/(3 db) {-(Sqrt[3]/2), -(1/2)};
k[3] = 4*Pi/(3 db) {(Sqrt[3]/2), -(1/2)};
a[1] = a[2] = a[3] = 1;


Now define the wave functions (this is where you had your problems - you had defined r as a constant, which isn't what your functions say)

ψs[r_] := Sum[a[i] Cos[(k[i]).r], {i, 3}];
ψa[r_] := Sum[a[i] Sin[(k[i]).r], {i, 3}];


And then it isn't clear whether you wanted to plot the square of the sum, $\left(\Psi_s(r)+\Psi_a(r)\right)^2$, or the sum of the squares, $\left(\Psi_s^2(r)+\Psi_a^2(r)^2\right)$. Only the square of their sum would show quantum interference effects. Here are both plots:

Grid[{{
DensityPlot[(ψs[{x, y}] + ψa[{x, y}])^2, {x, -10, 10}, {y, -10, 10}, PlotPoints -> 100, ImageSize -> 400],
DensityPlot[(ψs[{x, y}]^2 + ψa[{x, y}]^2), {x, -10, 10}, {y, -10, 10}, PlotPoints -> 100, ImageSize -> 400]
}}]


Edit: All of the tips given by Verbeia are very important, I didn't add them here so as not to repeat her - I just recognized this was in my area of expertise and added what I could.

• I see, this makes a lot more sense. When I defined r as {1,1,0}, aren't the cartesian unit vectors x,y,z already built in though when I use the curly brackets? Otherwise thank you, this helped immensely, I'm new with Mathematica so this was very enlightening :) – Rosie Sep 30 '15 at 20:51
• Thanks @JasonB - I should point out, though, that I'm a girl :-) The physics bit is not my area of expertise, so thank you for the complementary answer. – Verbeia Sep 30 '15 at 22:38

There are several things going on with your existing code. The first tip is that using Subscript can cause complications, so people tend to recommend that beginners don't use it for variable names.

The second tip is that the VectorAnalysis package is obsolete since version 9, and in any case, the usual Dot (.) function built into the main system will do what you want.

The third (and more fundamental) issue is that you have not actually defined functions of x and y for your DensityPlot to address. The x and y used in a plotting function are not the same as the ones you use outside that plot to define the expression to be plotted. This is a scoping issue; I recommend you have a look at some of the documentation on this issue. You instead need to define a function that takes a parameter, using the SetDelayed notation (:=). You would, for example, rewrite your first $\psi$ function as:

psis[x_] :=
Cos[DotProduct[k1, r]]*x + Cos[DotProduct[k2, r]]*x + Cos[DotProduct[k3, r]]*x;


Or using the built-in Dot:

psis1[x_] := Cos[k1.r]*x + Cos[k2.r]*x + Cos[k3.r]*x


Re-writing the function shows another issue: your functions are linear in x and y respectively. I think you actually want to rewrite them as follows.

psis2[x_] := Cos[k1.r*x] + Cos[k2.r*x] + Cos[k3.r*x];
psia2[y_] := Sin[k1.r*y] + Sin[k2.r*y] + Sin[k3.r*y];


This becomes obvious when you plot each version (note I've added something to psis1 because it's otherwise identical to psis):

Plot[{psis[x], psis1[x] + 10, 5 psis2[x]}, {x, 0, 100}]


So using these changed versions of your functions, we get closer to what you expect:

DensityPlot[(psis2[x])^2 + ( psia2[y])^2, {x, 0, 100}, {y, 0, 100},
ColorFunction -> GrayLevel]


It's not quite as in your desired picture, which leads me to suspect that you need to tweak your functional forms a bit more. But as I'm an economist and not a physicist, I'll have to leave that bit to you.

• Thank you, I will look into the resource you provided. What sort of problems can one encounter by using a subscript to define a variable? – Rosie Sep 30 '15 at 20:55
• Hi @Rosie - as noted here, if you see here, you will see that subscripted symbols are not plain symbols, so you can’t assign values (strictly speaking, DownValues`) to them directly – Verbeia Sep 30 '15 at 22:45