# Plotting Hydrogen's $2P_{x,y,z}$ Probability Densities in MATLAB [closed]

I have spent an unreasonable amount of time trying to plot $F(r,\theta,\phi)$ plane slices in MATLAB. I want to look at $x-y,y-z,x-z$ planes. Here's the function, specifically:

$$F(\rho,\theta,\phi)~=~4\pi\rho^{2}|\frac{\rho e^{\frac{-\rho}{2}}}{4\sqrt[]{2\pi}}\cos(\theta)|^{2}.$$

I certainly know about contour() and imagesc()in MATLAB. I would appreciate it if people would not respond "use contour()/imagesc()".

Unnecessary backstory: There happens to be a lot of information out there about plotting in spherical coordinates, however it's not relevant for my specific case because generally people have functions like $\rho=\theta^{2}$ to plot.

For the curious helper: This is Hydrogen's probability density function. The argument of the absolute value $|arg|$ is the $2P_{z}$ normalized wave function $\Psi(\rho,\theta,\phi)$. The $4\pi\rho^{2}$ accounts for the probability density conversion when working in spherical coordinates. People generally leave off the $4\pi$ portion since it is a constant.

-

## closed as off topic by David Z♦Apr 16 '13 at 7:12

Questions on Physics Stack Exchange are expected to relate to physics within the scope defined by the community. Consider editing the question or leaving comments for improvement if you believe the question can be reworded to fit within the scope. Read more about reopening questions here.If this question can be reworded to fit the rules in the help center, please edit the question.

What exactly are you asking here? If it's a question about MATLAB usage, it probably belongs on Super User (if that's the case I can migrate it for you). –  David Z Nov 9 '11 at 17:48
Well, I think it's a hybrid question that is related to MATLAB. I'm good with migrating to whichever community can best handle this. –  nick_name Nov 9 '11 at 17:54
OK, good to know... but I also meant that it's not clear (to me, at least) exactly what you're having problems with. Are you unsure about how to manipulate the function to get the numerical values of $F$ for a particular plane? If so, that's something we could handle here. But if you're unsure about what commands to use to display those numerical values once you have them, that belongs on Super User. –  David Z Nov 9 '11 at 18:23
I'm unsure about how to manipulate the function to get the numerical values of $F$ for a particular plane. –  nick_name Nov 9 '11 at 18:55
@DavidZaslavsky Also, I'm not sure what the standard $2P_{x}$ and $2P_{y}$ equations are for the Hydrogen atom, since $<2,1,1|$ and $<2,1,-1|$ are represented by the same wave function with a plus/minus phi term. Is it $sin(\phi)^{2}$ for $2P_{x}$ and $cos(\phi)^{2}$ for $2P_{y}$? –  nick_name Nov 9 '11 at 18:58

Just to make sure we're all on the same page, I'll take a step back before taking two steps forward. The normalized wavefunctions for the $n = 2$ energy level are usually written as \begin{gather} \psi_{2,1,0}(r, \theta, \phi) = \left(\frac{1}{32\pi a_0^3}\right)^{1/2} \frac{r}{a_0} \mathrm{e}^{-r/2a_0} \cos(\theta), \\ \psi_{2,1,\pm1}(r, \theta, \phi) = \mp \left(\frac{1}{64\pi a_0^3}\right)^{1/2} \frac{r}{a_0} \mathrm{e}^{-r/2a_0} \sin(\theta) \mathrm{e}^{\pm\mathrm{i}\phi}, \end{gather} where $a_0 = \hbar^2/m_\mathrm{e}e$ is the Bohr radius and the subscripts refer to quantum numbers $n$, $l$, and $m$, respectively. Sometimes these are written in terms of a dimensionless radius $\rho = r/a_0$. Sometimes $a_0$ is just explicitly set to $1$, though once this is done it's hard to figure out how to put dimensions back in, so I'll keep it there explicitly.

The $2P$ orbital is just the equally-weighted sum of these three parts: $$\psi_{2P} = \frac{1}{\sqrt{3}} \left(\psi_{2,1,0} + \psi_{2,1,1} + \psi_{2,1,-1}\right).$$ Here I put in the $\sqrt{3}$ in order to keep the probability normalized to $1$. If you are only interested in electron densities, you can take the square magnitude at this point, which I'll call $F$ to try to bring notations back in line. This turns out to be $$F(r, \theta, \phi) = \left\lvert \psi_{2P}(r, \theta, \phi) \right\rvert^2 = \frac{1}{96\pi a_0^5} r^2 \mathrm{e}^{-r/a_0} \left(1 - \sin^2(\theta) \cos(2\phi)\right),$$ as you are free to check.1

Now, to plot these, you are probably going to simply evaluate this function at Cartesian pixels.2 Therefore you want this as a function of $x$, $y$, and $z$. The usual definitions of spherical coordinates tell us \begin{align} r & = \sqrt{x^2+y^2+z^2}, \\ \sin^2(\theta) & = \frac{z^2}{x^2+y^2+z^2}, \\ \cos(2\phi) & = \frac{x^2-y^2}{x^2+y^2}, \end{align} where the last identity is best seen by using the identity $\cos(2\phi) = \cos^2(\phi) - \sin^2(\phi)$.

At this point, simply plug the formulas back into the definition of $F$, setting one of $x$, $y$, or $z$ to $0$ at a time. For example, the $x{-}y$ slice corresponds to $z = 0$, so $$F_{x{-}y}(x, y) = \frac{x^2+y^2}{96\pi a_0^5} \mathrm{e}^{-\sqrt{x^2+y^2}/a_0} \left(1 - \frac{x^2-y^2}{(x^2+y^2)^2}\right).$$

That's all there is to the physics - the rest is just throwing this into the language of your choice.3

1 Students should check this, to get a better handle on efficiently spotting orthogonality relations when integrating well-chosen eigenfunctions of Hermitian operators. Only once you are old as me can you Mathematica your way out of this exercise ;)

2 That's certainly the standard behavior of imagesc(). My Matlab experience doesn't go much further - there are other stackexchanges for how to best code this.

3 In Matlab, a quick way would be something like

[X,Y] = meshgrid(-3:0.01:3);
F = (X.^2 + Y.^2) / (96*pi) * exp(-sqrt(X.^2+Y.^2)) * (1 - (X.^2-Y.^2)/(X.^2+Y.^2).^2);
imagesc(F)

-