A smooth path to plot hydrogen atom via Monte Carlo method

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 ﬁltered 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.


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 [3][4][5][6][7].
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 1 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.
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: 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 : 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 ∈ R | −r ≤ y ≤ 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 ∈ [−r, r]. 6. Generate a random w ∈ [−r, r]. 7. If w ∈ [−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 = (2r) 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.
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 2 , 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 [19][20][21]. Although for low dimensions cases, the fixed numerical quadratures 3 provide better convergence rates, this MCM takes advantage because its accuracy is almost independent of dimension [19,20]. 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 [21][22][23] 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 [21][22][23]25].
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 [21][22][23]. 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.

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 [10,24,25]. 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 [24,25].
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 [24][25][26]. 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 [26]. 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 [24,25]. 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].

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.
In the previous equation, the Hamiltonian operatorĤ is associated with classical mechanics by an analogy to the combined operators of kinetic and potential energy, while the term ψ is the eigenfunction or wavefunction, and E is the eigenvalue or energy, we get both simultaneously when solving this equation [33,34].
For the termsĤ 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 ψ.
The explanation comes from the concept of particlewave duality. The particles, by their nature, are located in points, on the other hand, the wavefunction (as the name suggests) spreads through space [33,34], and consequently, there is no reason to represent electrons in orbits. Born made an interpretation that tells us that ψ ( r, t) is a probability amplitude, although it can often be complex 4 , |ψ ( r, t)| 2 is probability density and |ψ ( r, t)| 2 d 3 r is the probability dP ( r, t) to find the particle in a volume d 3 r, located between r and r + d r, at a time t [33,34].
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%.

Monte Carlo method to provide density plots
A common strategy involving MCMs is a von Neumann acceptance-rejection technique [35]. In his method, pseudorandom 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 5 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.

The one-dimensional infinite potential well
The Hamiltonian previously introduced (see section 4) is explicitly given in the following is the potential term, is the reduced Planck constant, m is the mass of the particle, and ∇ 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: where V (x) is given by: 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 [31,36,37]): where n is a positive integer (n = 1, 2, 3, ...), also known as a quantum number, and it comes from the application of the boundary conditions of this problem. 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.
When we calculate the density of probability, we obtain the following expression.
In a similar way to the example of the circle, we will use the function |ψ n (x)| 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. The typical result is available in Figure 2.

The two-dimensional infinite potential well
In the two-dimensional stationary version of the same problem, the equation (3) is given by: The potential V (x, y) in this case is: where L x and L y are the dimensions of the well, and the solution of equation (10) is (see [36,38]): 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: From this point on, the one-dimensional case algorithm is reproduced again for the function |ψ nxny (x, y)| 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 ≤ |ψ nxny (x, y)| 2 .
The typical result is shown in Figure 3.

The three-dimensional infinite potential well
In the three-dimensional stationary version of the potential well, the equation (3) is given by: − 2m ∇ 2 + V (x, y, z) ψ (x, y, z) = Eψ (x, y, z) (14) where the potential V (x, y, z) is: and the solution of the equation (14) is (see [34,39]): 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: Once again, we will change the algorithm in the onedimensional 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 ≤ |ψ nxnynz (x, y, z)| 2 . We show a typical result in Figure 4:

The hydrogen atom
In previous problems involving particle in boxes, the potential V ( 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 6 , 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 7 . We show this relation mathematically in the following equation (see [32]): 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, φ e θ are the azimuthal and polar angles respectively (see Figure 5). The previous equation has the following analytical solution (see [40]): where the radial part R(r) is: (20) and, The term n is the radial quantum number, which allows values n = {1, 2, 3, ...}. The index l corresponds to the azimuthal quantum number, which can assume the values l = {0, .., n − 1}. The constant a 0 is known as the Bohr radius, Z is the atomic number (Z := 1 for the hydrogen atom), and the term L 2l+1 n+l (ρ) is the associated Laguerre polynomial, expressed as follows: 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].
In the angular part, we use the following relation: (23) where the constant m is an integer, it represents the magnetic quantum number and can assume the values m = {−l, ..., +l}. The term P |m| l (cos(θ)) comes from the next two relations (24) and (25), respectively known as polynomials of Legendre associated and Legendre polynomials.
At this point, we have all the terms of the 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: Note that in the equation (23) there is an exponential term that carries the imaginary number i, making ψ nlm (r, θ, φ) non-measurable, but |Y l m (θ, φ)| 2 does not contain the imaginary number 8 and independent of φ. So, |ψ nlm (r, θ, φ)| 2 is also real and is a function only of the coordinates r and θ. Now we got |ψ nlm (r, θ, φ)| 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.
r = x 2 + y 2 + z 2 θ = arc cos z For the domain of equation (26), we know that θ ∈ [0, π], and r ∈ [0, +∞). 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 ∆x, ∆y, ∆z, centralized at the origin of the coordinate system. Then we use the relation (27) to evaluate r L = (1/2) ∆x 2 + ∆y 2 + ∆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, |ψ nlm (r, θ, φ)| 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 9 into a 2D regular grid, over the r and θ "axis", and for each node (grid point) we calculate |ψ nlm (r, θ, φ)| 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.   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 |ψ n (x)| 2 in any way.
The final comment should be made about the Legendre associated polynomials presented in the 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. 4. Set an amount p-points to plot. 5. Start an integer i-counter and reset it to zero. 6. Generate a random x ∈ [−∆x/2, ∆x/2], y ∈ [−∆y/2, ∆y/2] e z ∈ [−∆z/2, ∆z/2]. 7. Evaluate r and θ as the functions of x, y e z. 8. Generate a random w ∈ [0, M ]. 9. If w ≤ |ψ nlm (r, θ, φ)| 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 |ψ nlm (r, θ, φ)| 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 θ 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 ∆z = 0, and so on. We show some orbitals in Figures 9, 10 and 11.

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.