Solving the one dimensional Schrödinger equation of

the hydrogen and helium atom using the finite

elements method
MOGAING FEUSI Verdiane Laure
Under the supervision of Prof [Link]
Master2 MCN, University of Strasbourg
February 8, 2024

Abstract ically solved. These algebraic equations take into

account the appropriate boundary conditions for
Solving the Schrödinger equation analytically is each element, as well as the interactions between
difficult, so in this paper, we will present the fi- the electrons in the case of helium.
nite element method to solve the radial part of This paper describes a finite element calculation
the Schrödinger equation for the hydrogen and he- of the Hydrogen and Helium energy and radial
lium atoms. The finite element method allows us wave [Link] will also compare our results with
to approximate the unknown wave functions us- those from the literature and other numerical meth-
ing polynomial interpolation, specifically the La- ods.
grange polynomial in our case. The numerical re-
sults for the ground states of the hydrogen and he-
lium wave functions, as well as the respective eigen- 2 Schrödinger’s equation for
values of -0.4999 a.u which lies 0.02% above the Bi- hydrogen and helium atom
nesh, Mowlavi and [Link] value and -2.8568
a.u which lies 1.61% above the [Link] and J.S and appropriate units
[Link] results can be improved by increasing
the number of selected points. We can therefore The Schrödinger equation independent of time is
conclude that the finite element method can be very written as follows :
useful and effective in solving physical problems. ℏ2 2
− ∇ ψ + V (x)ψ = ϵψ (1)
2m

1 Introduction In the following, we will write the Schrödinger

equation for hydrogen and helium in spherical co-
The finite elements method (FEM) offers alterna- ordinates, and then define the units.
tive approach for studying few-body systems by
alllowing one to obtain a direct solution of the 2.1 Schrödinger’s equation for Hy-
Schrödinger equation for both bound and con- drogen
tinium states[1].
The FEM is an algorithm that uses local in- We have:
terpolation methods to solve second order differ-
ℏ2 2 e2
ential equations describing boundary-value prob- − ∇ ψ− ψ = ϵψ (2)
2m 4πϵ0
lems[1][2]. Using this approximation, we reformu-
late the Schrödinger equation into a system of lin- To start, let’s consider the Schrödinger equation for
ear algebraic equations, which can then be numer- the hydrogen atom in spherical coordinates. The

1
wave function of the electron in this coordinate sys- This leads us to the following two equations:
tem can be written as:
−1 ⃗2 ∥ψ2 (r⃗2 )∥2
Z
z
∇1 ψ1 (r⃗1 ) − ϕ1 (r⃗1 ) + dr⃗2 ψ1 r⃗1 =
ψ(r, θ, ϕ) = R(r)γ(θ, ϕ) (3) 2 r1 ∥r⃗1 − r⃗2 ∥
(6)
(2)
With ψ= the wave function; r= the radial distance (ϵ − ϵ a )ψ (
1 1r
⃗ ) = E ψ
1 1 1( r
⃗ )
from the nucleus; θ= the polar angle measured from −1 2
⃗2 ψ (r⃗ ) − ψ (r⃗ ) + dr⃗ ∥ψ1 (r⃗1 )∥ ψ r⃗ =
Z
z
the z-axis, and ϕ= the azimuthal angle measured ∇ 2 2 2 2 1 2 2
2 2 r2 ∥r⃗1 − r⃗2 ∥
from the x-axis. (7)
Taking these points into consideration, we obtain: (ϵ − ϵ(1) )ψ (
2 2r
⃗ ) = E ψ
2 1 2( r
⃗ )
a

1 d 2 dR 2mr2 Ze2 With

(r )+ (ϵ + )+ (4) Z
R dr dr ℏ2 4πϵ0
ϵa = d⃗ ⃗2 − Z ∥ψ (⃗
ri )∇
ri [ψi (⃗ 2
i i ri )∥ ]
ri
1 ∂ ∂γ 1 ∂γ 2
sin(θ) (sin(θ )) + 2 =0 We note:
γ ∂θ ∂θ γ sin (θ) ∂ϕ2
|ψ(r⃗′ )|2
Z
−34 −31 VH (⃗r) = dr⃗′ (8)
ℏ = 1.054.10 J.s; m = 9.1.10 Kg; |⃗r − r⃗′ |
e = 1.602.10−19 C; The Hartree potential represents the electrostatic
potential energy associated with the interaction be-
ϵ0 = 8.854.10−12 F/m tween an electron and an atomic nucleus, taking
into account the influence of other electrons in the
We will use atomic units in the rest of our work.
system. The total energy of helium atom is:
1 m e2 Z
= 2 , a0 = 0.529.10−9 m. ϵ = 2E − drVH (r)|U (r)|2 (9)
a0 ℏ 4πϵ0

ℏ2 -It is important to note that these two equa-

1Ha = = 27.20eV tions(6,7) are similar to those of the hydrogen
a20 m
atom. The most important term is just the pres-
The separation of the radial part from the angular ence of the Hartree energy, which couples the two
part and in atomic units leads to the following wave functions.
equations: -Given that an electronic orbital is occupied by 2
electrons (spin 1/2 and -1/2), in practice we only
1 d2 U l(l + 1) U (r)
+ U (r) − = ϵU (r) (5) need to solve a single equation, as the wave function
2 dr2 r2 r is the same for both electrons.
With U(r)= rR(r)
3 Finite elements method
2.2 Schrödinger’s equation for He-
lium we approximate the function ψ(x) by:
N
Helium, having two electrons, requires solving its X
equation using an approximation called the Hartree ψ(x) = αi ϕi (x) (10)
i=0
approximation, which assumes that each electron
evolves independently in a potential created by the On each interval (see figure 1 below), we have:
atomic nucleus and the other electron. This means  xi − xi−1
that the quantum state of each electron is deter-  if x ∈ [xi−1 , xi ]
xi − xi−1


mined by completely ignoring the effect of the other 

electron. Thus, we assume: ϕi (x) = xi+1 − x
if x ∈ [xi , xi+1 ]



 xi+1 − xi
ψ(r⃗1 , r⃗2 ) = ψ(r⃗1 )ψ(r⃗2 )

0 if x < xi−1 or x > xi−1

2
 have N-1 equations to solve. The two remaining
equations are found by studying the boundary con-
ditions.

3.2 FEM for Schrödinger’s equation

of helium
Figure 1: finite element discretization As we mentioned in section 2.2, the Schrödinger
equation for helium is similar to that of hydrogen.
The major problem that arises is the determination
of the Hatree potential, which we denote as VH (r).
We will now reduce the Schrödinger’s equations For the calculation of the Hartree potential, we
for hydrogen(6) and helium(7) using this approxi- use the Poisson’s equation in the following form in
mation. atomic units: ∇ ⃗2 VH (⃗r) = −4πρ(⃗r).
ρ(r) represents the volumetric charge density in
C/m3 (SI units). In considering the charge den-
3.1 FEM for Schrödinger’s equation sity inside a sphere, we find that:ρ(r) = |u(r)|22 .
4πr
of Hydrogen ultimately, VH is calculated by solving the equa-
By multiplying the entire equation by K(x), which tion:
d2 |u(r)|2
has the same expression as ψ(x) and by integrating (rV H (r)) = − . (14)
2
the entire equation(7) between x0 and xN we |dr {z } | {zr }
obtain: terme1 terme2

Z xN Z xN Term 1 is computed in finite element using the ex-

d2 ψ(x)
dxK(x) +2ϵ dxK(x)ψ(x) (11) pression of Tk (x), and term 2 is computed in finite
x0 dx2 x0 element using the expression of Vk (x).
| {z } | {z }
T(x) E(x) It is important to note that the electrons repel each
Z xN other, at very large distances, the interaction be-
−2 K(x)v(x)ψ(x) = 0 tween them decreases but does not become zero.
x0 so, the equation to be solved in this case is in the
| {z }
V(x) form:
with T|β>=|g> + |K>
l(l + 1) 1
v(x) = − (12)
x2 x 3.2.1 Matrix representation
T(x), E(x), V(x) are written in FEM as:
 αk−1 1 1 We obtain a system of N+1 equations of the form
 T k (x) = − ( + )α k (T+V)|α> =-2Eϵ|α> to solve in order to find the
xk − xk−1 xk − xk−1 xk+1 − xk


α coefficients and therefore the eigenvalues and



 αk+1
+ eigenvectors corresponding to the hydrogen atom.


xk+1 − xk


The corresponding matrices for each equation of



1 1 system (13) and for vectors |g>,|K> are:
E(x) = α k−1 (x k − x k−1 ) + α k (x k+1 − x k−1 )+


 6 3

 1

 αk+1 (xk+1 − xk )
6
  
1 0 0 0 ...




 1 0 T11 T12 . . . . . .
 Vk (x) = (xk+1 − xk−1 )v(xk )αk
2 
0 T21 T22 T23 . . .

(13) 
0 0 T32 T33 . . .

 
 .. .. .. .. .. 
with k ∈ [1, N − 1] In principle, we have N+1 un- . . . . . 
knowns α. We observe that, from this system, we 0 0 ... ... 1

3
 5 Numerical results for finite
 
1 0 0 0 ...
0 E11 E12 ... . . .

0 E21 E22 E23

. . . elements methods
 
0 0 E32 E33 . . .

 .. .. .... ..

 To materialize our numerical approach, I have cho-
. . .. .  sen to use the Python programming language.
0 0 ... ... 1 Given that the orbital radius of helium is larger
  than that of hydrogen, we have chosen the box size
1 0 0 0...
0 V11 and number of points for hydrogen to be xmax=20,
0 0. . .

0 0 V22
 N=600, and for helium to be xmax=40 and N=900.
0. . .

0 0
 As we have highlighted in section 3.2.2, due to the
0 V32 . . .

 .. .. .... ..

 divergence of the potential as the orbital radius
. . .. .  approaches zero, we have U(0)=U(+∞)=0, there-
0 0 ... ... 1 fore ψ(0)=ψ(+∞)=0. The energy value for l=0,
   which is the ground state of hydrogen, is -0.4999
0 β0
a.u, approximately -13.6030 eV, and for helium it
g0   −t 10 β0

 
g1  −tN −2,N −1 βN −1 
 is -2.8568 a.u, approximately -77.7380 eV. These re-
 
 ..   ..

 sults are in good agreement with analytical calcu-
. .  lations in quantum mechanics . The wave functions
1 βN −1 of hydrogen and helium in their ground states are
represented in Figure 2. Also, the interaction en-
3.2.2 Boundary conditions ergy between the electrons and the presence proba-
bility function in Coulomb potential wells for these
We have: U(r)=rR(r). Let’s study the behavior of atoms is plotted in Figure (a),(b) and (c).
the electron’s movement around the nucleus.
-when the electron is far from the
nucleus,U(+∞)=0, so ψ(+∞)=0 6 analysis and interpretation
-However, when it is very close to the nucleus, that
is, when r approaches 0, we need to develop R(r) -We notice that the curve of the wave function
as a powerP series. of helium is steeper than that of hydrogen. This
∞
R(r) = rα i=0 ai ri is due to the presence of two electrons in the 1S
we substitute it into equation (5), and we obtain orbital of the helium atom, resulting in a strong
in the end: attraction between the electron and the nucleus.
-The energy of the ground state of helium gives
 ∞
X us -2.8568 a.u, which corresponds to an error of
a0 + ai ri+1 if


1.6% of the result obtained by F.S. Levin and



 i=0
[Link][4]. The main reason for this is the


 l = 0, r− > 0, = 0

U (r) = rR(r) = Hartree approximation, which neglects Pauli’s
∞
 X exclusion principle.
a0 rl+1 + ai ri+1 if


-We also observe that the number of chosen points



 i=0
plays an important role in the accuracy of the



l = 0, r− > 0, = 0 results. The higher the number, the more precise
the results become.
consequently, ψ(0)=0
-The same applies to the chosen box size, which
Overall, we have: ψ(0)=ψ(+∞)=0.
plays a crucial role in visualizing the behavior of
the wave functions. A restricted size does not
accurately predict the electronic properties of
4 Problem restatement these atoms using the finite element method.
We have compiled in tables 1 and 2 some studies
The entire algorithm for solving our problem is on the solution of the Schrödinger equation using
rephrased in the appendix. various methods, as well as an estimation of the

4
 error compared to my results respectively for
hydrogen and helium. In these table,RK stands
for Runge Kutta Method, FDM stands for finite
difference method and HF:Hartree-Fock.

Author’s name energy(ev) error

FEM(our result) -13.6030
Binesh,Mowlavi,Arabshah[3](RK) - 13.6056 0.019%

Table 1: Comparison between my result and the

one obtained by the Runge-Kutta method for the
ground state of hydrogen.

(a) Results obtained from numerical Finite elements Author’s name energy(a.u) error
method for the ground state wave functions of hydro- Hawk-Hardcastle[5](FDM) 2,9036 1.6118%
gen and Helium [Link] and J. S[4]FEM -2,9032 1.5982%
Roothaan and OI.[6](HF) -2.8617 0.1712%

Table 2: Some results from other studies and com-

parison with the result found for helium.

7 Conclusion
These calculations for hydrogen and helium have
shown that the finite element method is a precise
and alternative approach for solving physical prob-
[Link] numerical results of this study are in good
agreement with the analytical results of quantum
mechanics. For more precise values, we need to
(b) interaction energy between the electrons of Helium increase the number of selected points. However,
this leads to longer computation time. Addition-
ally, this method requires significant memory.

8 Perspective
What will we obtain for an atom with more than 2
electrons, such as lithium?
It is evident that the problem will be more complex
and will require a more refined approximation, such
as the Hartree-Fock approximation, which takes
into account the electrostatic interactions between
electrons and considers each electron to evolve in
an average field created by all other electrons.

(c) presence probability function

5
 APPENDIX

A Problem restatement
Calculation of the 1s orbital of the hydrogen atom

calculation of the Hartree potential associated with the charge density of the 1s orbital of hydrogen

Solving the Schrödinger equation for helium (1S orbital)

Calculation of the total energy of the helium atom ϵold = ϵ, ϵ=2E- drVH (r)|u(r)|2
R

comparison of ϵ and ϵo ld, |ϵ-ϵold |<tolérance

d2 |u( r)|2 |
calculation of the Hartree potential dr 2 (rvH (r)) =− r

end

B Used code

6
7
8
References
[1] K. J. Bathe(1982)
Finite Element Procedures in Engineering Analysis .
[2] K. J. [Link] and E. Wilson(1976)
Numerical Methods in Finite Ele- ment Analysis.

[3] A. Binesh, A.A. Mowlavi and H. Arabshah(2010)

Application of Runge-Kutta numerical method to solve the Schrödinger equation for hy-
drogen and positronium atoms.
[4] I.L. Hawk and D.L. Hardcastle(1978)
Finite-difference solution to the Schrödinger equation for the Helium isoelectronique se-
quence.
[5] C. C.J. Roothaan,Lester [Link] and A. W. Weiss(1960)
Analytical self-Consistent fiels functions for the atomic configurations 1s2 , 1s2 2s, and
1s2 2s2 .