A smooth path to plot hydrogen atom via Monte Carlo method

Abstract

In this paper, we show how to build a basic computer program using the Monte Carlo method to display the hydrogen atomic orbitals. For this, in a heuristic way, we applied a von Neumann acceptance-rejection method in simple problems of potential wells, and we end with the hydrogen orbitals representation. In this technique, we spread points uniformly on the 1D and 2D charts of probability density distributions, then we filtered points under these “curves or surfaces”, and we extended this logic to 3D cases. Throughout the work, we also made some comments to help beginner students better understand the term “wave function” present in the Schrödinger equation. Also, we made all source code available at a third-party platform, for any purpose under the MIT license.

Keywords:
hydrogen atom; quantum mechanics; Monte Carlo method; Neumann acceptance-rejection method

1. Introduction

The introduction of quantum mechanics, even in physics or chemistry undergraduate courses is usually done by a historical review of several atomic models. Bohr developed one of these models, and it was essential to explain concepts of energy quantization and emission of spectral lines predicted by Rydberg ^[1]. Although the Bohr model is obsolete, it can still be used for beginner students to get acquainted with the quantization theme, and this is convenient because this model requires relatively simple conceptual and mathematical treatment as compared to its successor.

On the other hand, this practice may induce students to develop a misperception of the most modern and precise structures of quantum mechanics ^[2]. Therefore, researchers in the teaching area have suggested using computational tools to explain the Schrödinger models, making the current concepts of quantum mechanics more understandable [³– ⁷].

In this context, the present study shows how to use a computer for learning in this area. In more detail, this project aims to teach undergraduate students to build their programs using the Monte Carlo technique in order to picture the hydrogen orbitals in a heuristic¹ way.

This method is convenient because it induces the student to understand the real meaning of probability density introduced by Born, which derives from the wave function present in the time-independent version of the Schrödinger equation. So, we hope that students could make a suitable link between this function and the concept of particle-wave duality. In that way, we will give a brief presentation of the Monte Carlo method and quantum mechanics, applying these in some simple systems and ending up with the hydrogen atom.

2. The Monte Carlo method - MCM

In many science fields like physics [⁶, ⁷, ⁹– ¹²], biology [⁹, ¹⁰, ¹³– ¹⁷], engineering ^[9], and finance ^[18] many problems are solved today by MCMs. That is, by a broad class of methods that use random sampling, to obtain numerical results or provide the probability distributions.

This simple formulation suggests we need a more proper definition, and indeed, this could come in handy at this point. However, as these MCMs have quickly become popular and new techniques, ideas, and concepts continue to appear, then it becomes difficult to establish an overview.

For this reason, we will introduce the MCM by applying it to the following problem: take, as an example, the boundary of an irregular picture on a rectangular sheet of paper. One way of estimating the inner region of this figure is by throwing many darts randomly inside this sheet. The area of this figure is approximately the area of paper multiplied by the ratio of darts that hit the drawing and any other point.

Let us illustrate the process using the method to calculate the area of the circle inscribed in a square, and also sketch these regions with points. In this case, the equation of the circumference centered on the origin is given by:

(1) $$r^2=x^2+y^2$$

where r is the radius of the circle, x and y are cartesian coordinates. Then we will make y dependent on x, keeping the r constant, so we obtain the following expression with domain $D_f$ and range $R_f$:

(2) $$\begin{array}{l}y(x)=\sqrt{r^2-x^2}\\ D_f:\{x\in ℝ\mid -r\le x\le r\}\\ R_f:\{y\in ℝ\mid 0\le y\le r\}\end{array}$$

This function is represented graphically by an arc above the x-axis. The bottom part is −y(x). In this case, the range of −y(x) will be $R_f:\{y\in ℝ\mid -r\le y\le 0\}$. So, we can use a sequence of steps below to draw the circle and get its area.

1. Set a positive r.

2. Set an amount p-points to plot.

3. Set an integer h-counter and reset it to zero.

4. Set an integer i-counter and reset it to zero.

5. Generate a random $x\in [-r,r]$.

6. Generate a random $w\in [-r,r]$.

7. If $w\in [-y(x),y(x)]$, then x and w are stored in the first list and add +1 to h. Otherwise, they are stored in the second list.

8. Add +1 to i.

9. Repeat all steps from the 5th until i equals p.

10. Evaluate a circle area by $S={\left(2r\right)}^{2}h/p$.

11. Make a single chart with the two lists with different colors.

Note that to chart the whole process can be summed up in randomly picking a pair of coordinates $\{x,w\}$, and then representing the points in the external and internal limits in different colors. See the result in Figure 1.

Figure 1
A circle (radius r = 3) inscribed in a square (side l = 2r), with two thousand points homogeneously distributed, covering the whole area. The top and bottom are given by the functions y(x) and −y(x) respectively.

Concerning the calculation of the circle area, the reader may be wondering: why not calculate it by integration? Indeed, the exact solution is usually the best choice. Nevertheless, there are cases where integrals are challenging to manipulate to get a known analytical solution², or even are multidimensional.

In these cases, we must use numerical tools such as Trapezoidal and Simpson rules or even use the MCM to solve integrals [¹⁹–²¹]. Although for low dimensions cases, the fixed numerical quadratures³ provide better convergence rates, this MCM takes advantage because its accuracy is almost independent of dimension [¹⁹, ²⁰]. Also, this MCM avoids technical issues that appear in the other alternatives ^[19]. For these reasons and also for being easily implementable, the MCM is seen as a robust numerical tool.

However, be careful because there are at least two ways to calculate areas or volumes by MCMs. The technique we applied in our algorithmic is called von Neumann acceptance-rejection method (also known as hit-or-miss method), while that MCM which works as quadrature numeric is called sample-mean, also called crude Monte Carlo (see [²¹– ²³] for disambiguation). By the way, comparing these two options, only the von Neumann method can sketch areas and volumes, while getting these results avoiding the whole calculus formulation [²¹– ²³, ²⁵].

Furthermore, for the area of simple figures like circles, it is easy to show that these two types of MCMs, the more we increase the sample size, the closer we get to the analytical solution [²¹– ²³]. Be that as it may, as we shall see further on, rather than get numerical results, we will use the MCM to fill multidimensional spaces with points.

3. Random numbers and sampling

For these MCMs to work is essential that generators provide an unpredictable sequence of numbers, and identically given by some probability distribution [¹⁰,²⁴,²⁵]. In other words, good random number generators must produce uncorrelated numbers. When random numbers are produced homogeneously in some range, we describe the generator as a uniform random number generator [²⁴, ²⁵].

Also, when it comes to a programming environment, pseudo-random numbers are widely used. That is: they seem to be random, but they appear in deterministic sequences initiated by seeds [²⁴– ²⁶]. However, for different sets of values to be provided, i.e., seeming unpredictable, the actual time given by the CPU clock is used at least once to choose a seed [²⁶]. That way, every time the programs run, they seem to give unique results.

Since random numbers are present at the core of MCMs, then they have some basic requirements that generators must satisfy in order to be considered good [²⁴,²⁵]. Some tests suggest that those present in the software are unreliable, but fortunately, there are adequate packages ready to be used in more popular programming languages ^[27].

4. The wave equation

After the introduction of the Bohr model, de Broglie, motivated by the photoelectric effect, suggested that there could be a correlation between the electromagnetic equations with particles in movement (see the English version of ^[28] in ^[29]). Subsequently, Schrödinger, inspired by de Broglie hypothesis, or the concept we today know as particle-wave duality, formulated his homonymous equation below, which was later interpreted by Born, and which currently represents the primary tool of all non-relativistic quantum mechanics.

(3) $$\widehat{H}\psi =E\psi$$

In the previous equation, the Hamiltonian operator $\widehat{H}$ is associated with classical mechanics by an analogy to the combined operators of kinetic and potential energy, while the term $\psi$ is the eigenfunction or wavefunction, and E is the eigenvalue or energy, we get both simultaneously when solving this equation [³³, ³⁴].

For the terms $\widehat{H}$ and E, if beginners in quantum mechanics understand the concept of Newtonian mechanics and energy quantization introduced by Bohr, then they should have no difficulty in assimilating them. However, if the idea of electrons moving around the nucleus in precise orbits persists, there is no way to make a correct connection with $\psi$.

The explanation comes from the concept of particle-wave duality. The particles, by their nature, are located in points, on the other hand, the wavefunction (as the name suggests) spreads through space [³³, ³⁴] and consequently, there is no reason to represent electrons in orbits. Born made an interpretation that tells us that $\psi \left(\overrightarrow{r},t\right)$ is a probability amplitude, although it can often be complex⁴, $\psi {\left(\overrightarrow{r},t\right)}^2$ is probability density and ${\left|\psi \left(\overrightarrow{r},t\right)\right|}^2{d}^{3}r$ is the probability $dP\left(\overrightarrow{r},t\right)$ to find the particle in a volume $d^{3}r$, located between $\overrightarrow{r}$ and $\overrightarrow{r}+d\overrightarrow{r}$, at a time t [³³, ³⁴].

(4) $${\left|\psi \left(\overrightarrow{r},t\right)\right|}^2{d}^{3}r=dP\left(\overrightarrow{r},t\right)$$

If we integrate the previous equation into the whole space S, then the particles will inevitably be somewhere. Therefore, the probability of finding them under these conditions is 100%.

(5) $${\displaystyle {\int }_S}{\left|\psi \left(\overrightarrow{r},t\right)\right|}^2{d}^{3}r=1$$

5. Monte Carlo method to provide density plots

A common strategy involving MCMs is a von Neumann acceptance-rejection technique ^[35]. In his method, pseudo-random numbers are generated uniformly to produce other non-uniform distributions. At some stage of this algorithm, the screening of these numbers is made using relational operators, which result in the establishment of acceptance and rejection zones, similarly⁵ to the example of the circle already introduced. Also, as we propose in this paper, this same von Neumann technique was used before to draw hydrogen orbitals ^[6].

From now on, we will apply this MCM, to produce images of the qualitative probability density, and the number of points should be in a subjective range, so that image is not empty (without contrast) or overloaded. We will start from typical potential well problems, and end up with the hydrogen atom.

5.1. The one-dimensional infinite potential well

The Hamiltonian previously introduced (see section 4) is explicitly given in the following form $\widehat{H}=-\hslash {\nabla }^2/2m+V(\overrightarrow{r},t)$, where $V(\overrightarrow{r},t)$ is the potential term, $\hslash$ is the reduced Planck constant, m is the mass of the particle, and ${\nabla }^2$ is the Laplacian operator. Also, the time-independent version for equation (3) in one dimension, applied to the problem of the infinite potential well of length L is:

(6) $$-\frac{\hslash }{2m}\frac{d^2\psi \left(x\right)}{d{x}^2}+V\left(x\right)\psi \left(x\right)=E\psi \left(x\right)$$

where $V(x)$ is given by:

(7) $$V\left(x\right)=\{\begin{array}{ll}0,\hfill & \text{if}\text{0}<\text{x}<\text{L}\hfill \\ +\infty ,\hfill & \text{otherwise}.\hfill \end{array}$$

In this case, the particle can only occupy the limited region in the box from 0 to L, and the solution of the previous differential equation, including the application of the boundary conditions in this same interval is (see [³¹, ³⁶, ³⁷]):

(8) $${\psi }_n\left(x\right)=\sqrt{\frac{2}{L}}\sin \left(\frac{n\pi }{L}x\right)$$

where n is a positive integer $(n=1,2,3,\mathrm{...})$, also known as a quantum number, and it comes from the application of the boundary conditions of this problem. When we calculate the density of probability, we obtain the following expression.

(9) $$\begin{array}{l}{\left|{\psi }_n\left(x\right)\right|}^2=\frac{2}{L}{\sin }^2\left(\frac{n\pi }{L}x\right)\hfill \\ D_f:\left\{x\in ℝ\mid 0\le x\le L_x\right\}\hfill \\ R_f:\left\{{\left|{\psi }_n(x)\right|}^2\in ℝ\mid 0\le {\left|{\psi }_n(x)\right|}^2\le 2/L\right\}\hfill \end{array}$$

In a similar way to the example of the circle, we will use the function ${\left|{\psi }_n\left(x\right)\right|}^2$ as the filter to draw only the points within the domain, and ranges of the previous expression, by following the steps below:

1. Set a positive quantum number n.

2. Set an amount p-points to plot.

3. Set an integer i-counter and reset it to zero.

4. Generate a random $x\in [0,L]$.

5. Generate a random $w\in [0,2/L]$.

6. If $w\le {\left|{\psi }_n(x)\right|}^2$, then x and w are stored in the list and add +1 to i.

7. Repeat all steps from the 4th until i equals p.

8. Draw a graph with $\{x,w\}$ from the list.

9. Draw the curve ${\left|{\psi }_n(x)\right|}^2$ (optional).

The typical result is available in Figure 2.

Figure 2
Two thousand points below the probability density function ${\left|\psi (x)\right|}^2$, for a well of infinite potential of dimension L = 7, and quantum number n = 3.

5.2. The two-dimensional infinite potential well

In the two-dimensional stationary version of the same problem, the equation (3) is given by:

(10) $$\begin{array}{l}-\frac{\hslash }{2m}\left(\frac{{\partial }^2\psi \left(x,y\right)}{\partial x^2}+\frac{{\partial }^2\psi \left(x,y\right)}{\partial y^2}\right)V\left(x,y\right)\psi \left(x,y\right)\\ =E\psi \left(x,y\right)\end{array}$$

The potential $V(x,y)$ in this case is:

(11) $$V\left(x,y\right)=\{\begin{array}{ll}0,\hfill & \text{if}\text{0}\le \text{x}\le {\text{L}}_{\text{x}},\text{0}\le \text{y}\le {\text{L}}_{\text{y}}\hfill \\ +\infty ,\hfill & \text{otherwise}.\hfill \end{array}$$

where $L_x$ and $L_y$ are the dimensions of the well, and the solution of equation (10) is (see [³⁶, ³⁸]):

(12) $${\psi }_{n_x{n}_y}\left(x,y\right)=\frac{2}{\sqrt{L_x{L}_y}}\sin \left(\frac{n_x\pi }{L_x}x\right)\sin \left(\frac{n_y\pi }{L_y}y\right)$$

In the last expression, the numbers $n_x$ and $n_y$ are the positive integers that appear when the boundary conditions are applied. When calculating the probability density, we obtain the following relation:

(13) $$\begin{array}{l}{\left|{\psi }_{n_x{n}_y}\left(x,y\right)\right|}^2=\frac{4}{L_x{L}_y}{\sin }^2\left(\frac{n_x\pi }{L_x}x\right){\sin }^2\left(\frac{n_y\pi }{L_y}y\right)\\ D_f:\left\{(x,y)\in ℝ\mid 0\le x\le L_x,0\le y\le L_y\right\}\\ R_f:\left\{{\left|{\psi }_{n_x{n}_y}(x,y)\right|}^2\in ℝ\mid 0\le {\left|{\psi }_{n_x{n}_y}(x,y)\right|}^2\le \frac{4}{L_x{L}_y}\right\}\end{array}$$

From this point on, the one-dimensional case algorithm is reproduced again for the function ${\left|{\psi }_{n_x{n}_y}\left(x,y\right)\right|}^2$, making some adaptations. In the 1st item we set the values of $n_x$ and $n_y$, in the 4th and 5th items, we generate the x, y, and w in the intervals defined in domain $D_f$ and range $R_f$, and then we make changes in the 6th and 8th items to plot the points that satisfy $w\le {\left|{\psi }_{n_x{n}_y}\left(x,y\right)\right|}^2$. The typical result is shown in Figure 3.

Figure 3
Five thousand points distributed below the probability density function ${\left|\psi (x,y)\right|}^2$, for a infinite potential well of dimensions $L_x=25$, $L_y=30$, and quantum numbers $n_x=3$, and $n_y=2$.

5.3. The three-dimensional infinite potential well

In the three-dimensional stationary version of the potential well, the equation (3) is given by:

(14) $$\left(-\frac{\hslash }{2m}{\nabla }^2+V\left(x,y,z\right)\right)\psi \left(x,y,z\right)=E\psi \left(x,y,z\right)$$

where the potential $V(x,y,z)$ is:

(15) $$V\left(x,y,z\right)=\{\begin{array}{ll}0,\hfill & \text{if}\text{0}\le \text{x}\le {\text{L}}_{\text{x}},\hfill \\ \hfill & 0\le y\le L_y,0\le z\le L_z\hfill \\ +\infty ,\hfill & \text{otherwise}.\hfill \end{array}$$

and the solution of the equation (14) is (see [³⁴, ³⁹]:

(16) $$\begin{array}{lll}{\psi }_{n_x{n}_y{n}_z}\left(x,y,z\right)\hfill & =\hfill & \sqrt{\frac{8}{L_x{L}_y{L}_z}}\sin \left(\frac{n_x\pi }{L_x}x\right)\hfill \\ \hfill & \hfill & \sin \left(\frac{n_y\pi }{L_y}y\right)\sin \left(\frac{n_z\pi }{L_z}z\right)\hfill \end{array}$$

The constant $L_z$ is the additional dimension of the box, and $n_z$ is the third quantum number. When we calculate the density of probability, we obtain the following expression:

(17) $$\begin{array}{l}\mid {\psi }_{n_x{n}_y{n}_z}(x,y,z){\mid }^2=\frac{8}{L_x{L}_y{L}_z}\hfill \\ {\sin }^2(\frac{n_x\pi }{L_x}x){\sin }^2(\frac{n_y\pi }{L_y}y){\sin }^2(\frac{n_z\pi }{L_z}z)\hfill \\ D_f:\{(x,y,z)\in ℝ\mid 0\le x\le L_x,0\le y\le L_y,0\le z\le L_z\}\hfill \\ R_f:\{\mid {\psi }_{n_x{n}_y{n}_z}(x,y,z){\mid }^2\in ℝ\mid 0\hfill \\ \le \mid {\psi }_{n_x{n}_y{n}_z}(x,y,z){\mid }^2\le \frac{8}{L_x{L}_y{L}_z}\}\hfill \end{array}$$

Once again, we will change the algorithm in the one-dimensional case to include new coordinates, box dimensions, and quantum numbers. In this case, the program will provide x, y, z, and the term w. Although, there is no way to graph these four variables simultaneously, and we will represent the points that satisfy $w\le {\left|{\psi }_{n_x{n}_y{n}_z}(x,y,z)\right|}^2$. We show a typical result in Figure 4:

Figure 4
Fifty thousand points distributed in three-dimensional space, whose values of w are below the probability density function ${\left|\psi (x,y,z)\right|}^2$, for a infinite potential well, with dimensions $L_x=90,L_y=30,L_z=60$, and quantum numbers $n_x=2,n_y=1,n_z=4$.

5.4. The hydrogen atom

In previous problems involving particle in boxes, the potential $V(\overrightarrow{r})$ is responsible for the confinement of the particle, being infinite outside of limits and null in its interior. However, for the hydrogen atom alone⁶, i.e., a system with a single proton and its electron, the Coulomb force between them combined with the kinetic energy is responsible for the entire motion of the electron around the nucleus⁷. We show this relation mathematically in the following equation (see ^[32]):

(18) $$-\left(\frac{\hslash }{2m}{\nabla }^2+\frac{e^2}{4\pi {ϵ}_{0}r}\right)\psi \left(r,\theta ,\varphi \right)=E\psi \left(r,\theta ,\varphi \right)$$

The terms e and m are charge and mass of the electron, ${ϵ}_0$ is the constant of electrical permittivity in the vacuum. In this problem, we use spherical coordinates, i.e., r is the radius, $\varphi$ e $\theta$ are the azimuthal and polar angles respectively (see Figure 5). The previous equation has the following analytical solution (see ^[40]):

Figure 5
Spherical coordinates, where the radial distance is r, the polar angle is $\theta$, and the azimuth angle is $\varphi$.

(19) $$\psi \left(r,\theta ,\varphi \right)=R(r)Y(\theta ,\varphi )$$

where the radial part $R(r)$ is:

(20) $$R_{nl}(r)=-\sqrt{{\left(\frac{2Z}{n{a}_0}\right)}^3\frac{\left(n-l-1\right)!}{2n{\left[\left(n+l\right)!\right]}^3}}{e}^{-\rho /2}{\rho }^l{L}_{n+l}^{2l+1}\left(\rho \right)$$

and,

(21) $$\rho (r)\text{:=}\frac{2Zr}{n{a}_0}.$$

The term n is the radial quantum number, which allows values $n=\{1,2,3,\mathrm{...}\}$. The index l corresponds to the azimuthal quantum number, which can assume the values $l=\{0,\mathrm{..},n-1\}$. The constant $a_0$ is known as the Bohr radius, Z is the atomic number ($Z\text{:=}1$ for the hydrogen atom), and the term $L_{n+l}^{2l+1}(\rho )$ is the associated Laguerre polynomial, expressed as follows:

(22) $$\begin{array}{lll}{L}_{n+l}^{2l+1}\left(\rho \right)\hfill & =\hfill & {\displaystyle \sum _{k=0}^{n-l-1}{\left(-1\right)}^{k+1}}\hfill \\ \hfill & \hfill & \frac{{\left[\left(n+l\right)!\right]}^2}{\left(n-l-1-k\right)!\left(2l+1+k\right)!k!}{\rho }^k\hfill \end{array}$$

In the angular part, we use the following relation:

(23) $$Y_l^m(\theta ,\varphi )=\sqrt{\frac{\left(2l+1\right)\left(l-\left|m\right|\right)!}{4\pi \left(l+\left|m\right|\right)!}}{P}_l^{\left|m\right|}\left(\cos \left(\theta \right)\right)e^{im\varphi }$$

where the constant m is an integer, it represents the magnetic quantum number and can assume the values $m=\{-l,\mathrm{...},+l\}$. The term $P_l^{\left|m\right|}(\cos (\theta ))$ comes from the next two relations (24) and (25), respectively known as polynomials of Legendre associated and Legendre polynomials.

(24) $$P_l^{\left|m\right|}(x)={\left(-1\right)}^{\left|m\right|}{\left(1-x^2\right)}^{\left|m\right|/2}\frac{d^{\left|m\right|}}{d{x}^{\left|m\right|}}{P}_l(x)$$

(25) $$P_l(x)=\frac{1}{2^{l}l!}\frac{d^l}{d{x}^l}{\left(x^2-1\right)}^l$$

At this point, we have all the terms of the equation equation (19) available, and we can proceed by calculating the probability density numerically. Moreover, as an alternative, at least for the first quantum numbers, it is also possible to perform the calculations by the wave functions present in Table 1. Therefore, whichever path we choose, we must obtain the following probability density:

Table 1
First normalized wave functions of the hydrogen atom. This is an adaptation of the table 6.1 from ^[31].

n	l	m	${\psi }_{n,l,m}\left(r,\theta ,\varphi \right)$
1	0	0	$\frac{1}{\sqrt{\pi }{a}_0^{3/2}}{e}^{-r/a_0}$
2	0	0	$\frac{1}{4\sqrt{2\pi }{a}_0^{3/2}}\left(2-\frac{r}{a_0}\right)e^{-r/2a_0}$
2	1	0	$\frac{1}{4\sqrt{2\pi }{a}_0^{3/2}}\frac{r}{a_0}{e}^{-r/2a_0}cos\theta$
2	1	$\pm 1$	$\frac{1}{8\sqrt{3\pi }{a}_0^{3/2}}\frac{r}{a_0}{e}^{-r/2a_0}sin\theta e^{\pm i\varphi }$
3	0	0	$\frac{1}{81\sqrt{3\pi }{a}_0^{3/2}}\left(27-18\frac{r}{a_0}+2\frac{r^2}{a_o^2}\right)e^{-r/2a_o}$
3	1	0	$\frac{1}{81\sqrt{3\pi }{a}_0^{3/2}}\left(6-\frac{r}{a_0}\right)\frac{r}{a_0}{e}^{-r/3a_0}cos\theta$
3	1	$\pm 1$	$\frac{1}{81\sqrt{3\pi }{a}_0^{3/2}}\left(6-\frac{r}{a_0}\right)\frac{r}{a_0}{e}^{-r/3a_0}sin\theta e^{\pm i\varphi }$
3	2	0	$\frac{1}{81\sqrt{6\pi }{a}_0^{3/2}}\frac{r^2}{a_0^2}{e}^{-r/3a_0}\left(3co{s}^2\theta -1\right)$
3	2	$\pm 1$	$\frac{1}{81\sqrt{\pi }{a}_0^{3/2}}\frac{r^2}{a_0^2}{e}^{-r/3a_0}sin\theta cos\theta e^{\pm i\varphi }$
3	2	$\pm 2$	$\frac{1}{162\sqrt{\pi }{a}_0^{3/2}}\frac{r^2}{a_0^2}{e}^{-r/3a_0}si{n}^2\theta e^{\pm 2i\varphi }$

(26) $${\left|{\psi }_{nlm}\left(r,\theta ,\varphi \right)\right|}^2={\left|R_{nl}(r)\right|}^2{\left|Y_m^l(\theta ,\varphi )\right|}^2$$

Note that in the equation equation (23) there is an exponential term that carries the imaginary number i, making ${\psi }_{nlm}\left(r,\theta ,\varphi \right)$ non-measurable, but ${\left|Y_m^l(\theta ,\varphi )\right|}^2$ does not contain the imaginary number⁸ and independent of $\varphi$. So, ${\left|{\psi }_{nlm}\left(r,\theta ,\varphi \right)\right|}^2$ is also real and is a function only of the coordinates r and $\theta$.

Now we got ${\left|{\psi }_{nlm}\left(r,\theta ,\varphi \right)\right|}^2$ already explained, and we can determine its domain and the range to apply the MCM, as we have done before. So, the first thing we need to do is to be able to switch from cartesian to a spherical coordinate system. To achieve this, we can extract this information from Figure 5 to get the relation (27) below.

(27) $$\begin{array}{l}r=\sqrt{x^2+y^2+z^2}\hfill \\ \theta =arccos\left(\frac{z}{\sqrt{x^2+y^2+z^2}}\right)\hfill \end{array}$$

For the domain of equation (26), we know that $\theta \in [0,\pi ]$, and $r\in [0,+\infty )$. Since the radius r has no upper limit, it is impossible to search the electron in the entire space. Alternatively, we will check all volume by an arbitrary box which contains the location most likely to find the single electron.

Therefore, we place a box of dimensions $\Delta x,\Delta y,\Delta z$, centralized at the origin of the coordinate system. Then we use the relation equation (27) to evaluate $r_L=(1/2)\sqrt{\Delta x^2+\Delta y^2+\Delta z^2}$. This $r_L$ is the magnitude of a vector that goes from the origin to any corner and will be used as the upper limit of the domain.

In respect of range, ${\left|{\psi }_{nlm}\left(r,\theta ,\varphi \right)\right|}^2$ is lower bounded by 0, and it reaches different top values depending on the quantum numbers. Unfortunately, finding a generic function that calculates the global maximum for each set involving all of them is no easy task. A first alternative would be to find a solution analytically for each case, but this process becomes increasingly tedious as quantum numbers grow. Another option would be to use some searching process to seek global maximum numerically ^[41], after all, we have the computer available. Instead, let us proceed more smoothly.

The simplest way is to split the space⁹ into a 2D regular grid, over the r and $\theta$ "axis", and for each node (grid point) we calculate ${\left|{\psi }_{nlm}\left(r,\theta ,\varphi \right)\right|}^2$. Then we update a variable that stores the highest value M found, and then increase it to ensure it is above the global maximum of this function. The Figures 6, 7, and 8 exemplify this method.

Figure 6
The probability density as a function of r and $\theta$ for the hydrogen atom, with quantum numbers $n=1,l=0,m=0$. The points are projections of nodes in the function ${\left|\psi (r,\theta ,\varphi )\right|}^2$.

Figure 7
The probability density as a function of r and $\theta$ for the hydrogen atom, with quantum numbers $n=3,l=2,m=0$. The points are projections of nodes in the function ${\left|\psi (r,\theta ,\varphi )\right|}^2$.

Figure 8
The probability density as a function of r and $\theta$ for the hydrogen atom, with quantum numbers $n=3,l=2,m=1$. The points are projections of nodes in the function ${\left|\psi (r,\theta ,\varphi )\right|}^2$.

This procedure may not be elegant, but it is practical. If the reader is uncomfortable about this, do the same with the problem in $1D$, then he will realize that the method will discard points above the function ${\left|{\psi }_n\left(x\right)\right|}^2$ in any way.

The final comment should be made about the Legendre associated polynomials presented in the equation equation (24). Some environments operate with algebraic or symbolic computations to interpret equations like that, but to reproduce this method in several programming languages, we must soften this step. In this case, it is better to use suitable libraries for this purpose, because the implementation of simpler algorithms may cause problems of numerical instability ^[42]. A good choice for the C++ programming language is the BOOST library ^[43] and the SHTOOLS^[44] for FORTRAN and Python.

Finally, with all this information gathered, we have the conditions to represent the orbitals of the hydrogen atom. For that, we will adapt the algorithm seen in the examples of the confined particle in boxes, presented in the following.

1. Set the quantum numbers $n\in [1,+\infty )$, $l\in [0,n-1]$, $m\in [-l,l]$.

2. Set the box dimensions $\Delta x$, $\Delta y$ e $\Delta z$.

3. Estimate an upper limit for M above of ${\left|{\psi }_{nlm}\left(r,\theta ,\varphi \right)\right|}^2$.

4. Set an amount p-points to plot.

5. Start an integer i-counter and reset it to zero.

6. Generate a random $x\in [-\Delta x/2,\Delta x/2],y\in [-\Delta y/2,\Delta y/2]$ e $z\in [-\Delta z/2,\Delta z/2]$.

7. Evaluate r and $\theta$ as the functions of x, y e z.

8. Generate a random $w\in [0,M]$.

9. If $w\le {\left|{\psi }_{nlm}(r,\theta ,\varphi )\right|}^2$, them x, y and z are stored in the list, and add +1 to i.

10. Repeat all steps from the 6th until i equals p.

11. Make a graph with the coordinates from the list.

To estimate the upper limit M, we calculate $r_L$ to set up the grid, then we scan it to get an approximate global maximum of ${\left|{\psi }_{nlm}\left(r,\theta ,\varphi \right)\right|}^2$, and finally, we add 5% over this value found. To make this grid, we suggest that use nearly 50 nodes per Angström along of the r radius, and 160 nodes over the $\theta$ angle.

For plotting the orbitals in 2D, turn the box into a rectangle by choosing the null value for some side. For example, if the interest is in the $x-y$ plane, then $\Delta z=0$, and so on. We show some orbitals in Figures 9, 10 and 11.

Figure 9
The simulated hydrogen atom with fifteen thousand points.

Figure 10
The simulated hydrogen atom with fifteen thousand points.

Figure 11
The simulated hydrogen atom with fifteen thousand points.

6. Conclusion

This paper shows in a heuristic way that it is easy to represent the orbitals of the hydrogen atom using basic quantum mechanics and coding skills.

We also made all source code available at a third-party platform ^[45], for any uses under the MIT license.

Acknowledgment

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

• 1
  Heuristic is “a method which, on the basis of experience or judgement, seems likely to yield a reasonable solution to a problem, but which cannot be guaranteed to produce the mathematically optimal solution.” – Edward A. Silver (see ^[8]).
• 2
  It is not always easy to substitute non-elementary integrals for simple power series anti-derivative functions.
• 3
  Numerical quadrature is nothing more than a technique for making approximations of integrals by sums.
• 4
  Complex in this context because it contains the imaginary number i.
• 5
  In his original work ^[35], von Neumann created a linear function to convert one interval of random numbers into another. We thought this was necessary because, the generators at the time, perhaps only would provide values from 0 to 1. Today there are tools in programming languages that provide random numbers in any desired range, and we can take advantage of this in some cases.
• 6
  The same applies to hydrogen-like ions such as ${\text{He}}^{+},\text{L}{\text{i}}^{2+},{\text{Be}}^{\text{3+}},{\text{B}}^{4+}$, and so on. The only difference is that the atomic number Z must be changed in each case.
• 7
  If we were strict, we should consider these two bodies could vibrate, rotate around the center mass of the system, and along to their axis, or even interact with externals fields. Although, all these effects are minimal and are usually neglected to ensure that the problem has an analytical solution ^[32].
• 8
  The product between a number and its conjugate always gives a non-negative real number: $z\cdot z^{*}={\left|z\right|}^2$
• 9
  For example, r and $\theta$ starts at 0 and reaches $r_L$ and $\pi$ respectively, so if the length $r_L$ is divided by 5 and angle $\pi$ by 4, then we have 30 nodes $(30=(5+1)\times (4+1))$ to check.

References

• ^[1]
  A. Lakhtakia, Models and Modelers of Hydrogen (World Scientific Publishing Company, Singapura 1996).
• ^[2]
  H. Fischler and M. Lichtfeldt, International Journal of Science Education 14, 181 (1992).
• ^[3]
  D.T. Cromer, J. Chem. Educ. 45, 626 (1968).
• ^[4]
  B. Coto, A. Arencibia and I. Suárez, Computer Applications in Engineering Education 24, 765 (2016).
• ^[5]
  B.G. Moore, J. Chem. Educ. 77, 785 (2000).
• ^[6]
  M. Goto and V.M. Aquino, Semina: Ciências Exatas e Tecnológicas 13, 255 (1992).
• ^[7]
  V.M. Aquino, V.C. Aguilera-Navarro, M. Goto and H. Iwamoto, American Journal of Physics 69, 788 (2001).
• ^[8]
  E.A. Silver, Journal of the Operational Research Society 55, 936 (2004).
• ^[9]
  M. Shlomo and S. Mordechai, Applications of Monte Carlo Method in Science and Engineering (IntechOpen, London, 2011).
• ^[10]
  D.P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, Second Edition (Cambridge University Press, New York, 2005).
• ^[11]
  V.L. Líbero, Rev. Bras. Ens. Fis. 22, 346 (2000).
• ^[12]
  R. da Silva and J.R.D. de Felício, Rev. Bras. Ens. Fis. 24, 103 (2002).
• ^[13]
  E. Arashiro, J.R.D. de Felício and U.H.E. Hansmann, Physical Review E 73, 040902 (2006).
• ^[14]
  F.M. Ruziska, E. Arashiro and T. Tomé, Physica A: Statistical Mechanics and its Applications 489, 56 (2018).
• ^[15]
  N.R.S. Ortega, C.F.S. Pinheiro, T. Tomé and J.R.D. de Felício, Physica A: Statistical Mechanics and its Applications 255, 189 (1998).
• ^[16]
  F. Keesen, A. Castro e Silva, E. Arashiro, C.F.S. Pinheiro, Ecological Modelling 344, 38 (2017).
• ^[17]
  A. Castro e Silva and A.T. Bernardes, Physica A: Statistical Mechanics and its Applications 352, 535 (2005).
• ^[18]
  C.A. Waldspurger, T. Hogg, B.A. Huberman, J.O. Kephart and W.S. Stornetta, IEEE Transactions on Software Engineering 18, 103 (1992).
• ^[19]
  F. James, Reports on Progress in Physics 43, 1145 (1980).
• ^[20]
  R.E. Caflisch, Acta Numerica 7, 1 (1998).
• ^[21]
  H. Gould, J. Tobochnik and W. Christian An introduction to computer simulation methods: applications to physical systems-3rd Edition (Pearson Addison Wesley, 2007).
• ^[22]
  R.Y. Rubinstein, Simulation and the Monte Carlo Method (Wiley, Hoboken, 1981).
• ^[23]
  Y. Badis, Computational Physics (WORLD SCIENTIFIC, Singapore, 2017).
• ^[24]
  D.P. Kroese, T. Taimre and Z.I. Botev, Handbook of Monte Carlo Methods (John Wiley & Sons, Hoboken, 2011).
• ^[25]
  J.K.S. William and L. Dunn, Exploring Monte Carlo Methods (Elsevier, Cambridge, 2011).
• ^[26]
  M. Lee, C++ Programming for the Absolute Beginner, 2nd Edition (Course Technology Press, Boston, 2009).
• ^[27]
  https://ieeexplore.ieee.org/document/977250
  » https://ieeexplore.ieee.org/document/977250
• ^[28]
  https://www.academie-sciences.fr/pdf/dossiers/Broglie/Broglie_pdf/CR1923_p507.pdf
  » https://www.academie-sciences.fr/pdf/dossiers/Broglie/Broglie_pdf/CR1923_p507.pdf
• ^[29]
  https://www.nature.com/articles/112540a0
  » https://www.nature.com/articles/112540a0
• ^[30]
  V.V. Raman and P. Forman, Historical Studies in the Physical Sciences 1, 291 (1969).
• ^[31]
  A. Beiser, Concepts of Modern Physics (McGraw-Hill, Boston, 2002).
• ^[32]
  I. Levine, The Hydrogen Atom (Pearson, Boston, 2014).
• ^[33]
  D.J. Griffiths Introduction to Quantum Mechanics (Cambridge University Press, Cambridge, 2017).
• ^[34]
  N. Zettili, Quantum Mechanics: Concepts and Applications (Wiley, Hoboken, 2009).
• ^[35]
  J. Von Neumann, National Bureau of Standards Applied Mathematics 3, 36 (1951).
• ^[36]
  G.D.W.W. Bauer, in University Physics with Modern Physics (McGraw-Hill, New York, 2011).
• ^[37]
  https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Quantum_Mechanics/05.5\%3A_Particle_in_Boxes/Particle_in_a_1-Dimensional_box
  » https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Quantum_Mechanics/05.5\%3A_Particle_in_Boxes/Particle_in_a_1-Dimensional_box
• ^[38]
  https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Quantum_Mechanics/05.5\%3A_Particle_in_Boxes/Particle_in_a_2-Dimensional_box
  » https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Quantum_Mechanics/05.5\%3A_Particle_in_Boxes/Particle_in_a_2-Dimensional_box
• ^[39]
  https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Quantum_Mechanics/05.5\%3A_Particle_in_Boxes/Particle_in_a_3-Dimensional_box
  » https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Quantum_Mechanics/05.5\%3A_Particle_in_Boxes/Particle_in_a_3-Dimensional_box
• ^[40]
  L. Pauling, Introduction to quantum mechanics: with applications to chemistry (Dover Publications, New York, 1985).
• ^[41]
  J. Nocedal and S. Wright Numerical Optimization (Springer Science & Business Media, Berlin, 2006).
• ^[42]
  W. Press, Numerical recipes in C : the art of scientific computing (Cambridge University Press, Cambridge, 1992).
• ^[43]
  https://www.boost.org/doc/libs/1_56_0/libs/math/doc/html/index.html
  » https://www.boost.org/doc/libs/1_56_0/libs/math/doc/html/index.html
• ^[44]
  https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018GC007529
  » https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018GC007529
• ^[45]
  https://github.com/pedrohflobo/hydrogen_orbitals_by_Monte_Carlo_Method
  » https://github.com/pedrohflobo/hydrogen_orbitals_by_Monte_Carlo_Method

Publication Dates

• Publication in this collection
  22 July 2019
• Date of issue
  2019

History

• Received
  30 Mar 2019
• Reviewed
  16 May 2019
• Accepted
  25 May 2019