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

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

enter image description here

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

enter image description here

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

  • 4
    $\begingroup$ 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. $\endgroup$
    – Verbeia
    Commented Sep 30, 2015 at 3:53
  • $\begingroup$ @Rosie, do not forget to tick (not only upvote) the answer if it solves your problem =)) $\endgroup$
    – garej
    Commented Sep 30, 2015 at 20:24
  • $\begingroup$ @garej Thanks I will do so $\endgroup$
    – Rosie
    Commented Sep 30, 2015 at 21:06

2 Answers 2


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:

  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]

density plots

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.

  • $\begingroup$ 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 :) $\endgroup$
    – Rosie
    Commented Sep 30, 2015 at 20:51
  • $\begingroup$ 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. $\endgroup$
    – Verbeia
    Commented Sep 30, 2015 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}]

enter image description here

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]

enter image description here

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.

  • $\begingroup$ 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? $\endgroup$
    – Rosie
    Commented Sep 30, 2015 at 20:55
  • $\begingroup$ 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 $\endgroup$
    – Verbeia
    Commented Sep 30, 2015 at 22:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.