Numerical solution of the Schrödinger equation using Neural Networks in Python

The motion of quantum mechanical systems in physical sciences is described by partial differential equations, usually of second order with respect to spatial coordinates. The required solutions of the time-dependent type of equations are, in general, functions of the temporal variable t and the state vectors determining the positions of the system’s particles at the time t. Only a small number, however of those differential equations can however be solved analytically, while the majority of them must be solved numerically by applying specific very advanced integration techniques. Among these equations, the fundamental Schrödinger equation offers great insight towards numerical solving. In this work we present an effective method that solves numerically the time-independent Schrodinger equation on the basis of neural networks techniques. Analytical and numerical results for the radial part of this equation with the corresponding energies are then compared so as to estimate the performance of our method.


Introduction
It is well known that, many phenomena in almost every scientific field including classical and quantum physics, can be described through the use of partial differential equations (PDEs), in the vast majority of which, finding an analytical solution is impossible.Many techniques, computational, numerical and heuristic have been developed in various attempts to model physical phenomena [1,2,3,4,5,6,7].Towards this aim, about a century ago, Schrödinger derived a second order differential equation to account for the wave nature of a particle assuming that the particle behaves as wave.This equation is known as a wave-equation and serves as a mathematical model to describe the motion of a particle by means of a wave [8,9].The Schrödinger's equation yields solutions in terms of functions, called wave functions, while at the same time the corresponding energy E of the particle in question is obtained.In general, the wave-functions can describe the motion of microscopic systems in molecular, atomic and subatomic (nuclear and hadronic) physics but also of macroscopic systems, for example, the entire Universe.Hence, Schrödinger's equation is considered fundamental to all applications of quantum mechanics including the quantum field theory.
In mathematics, Schrödinger's equation in its two variants: the time-dependent and the timeindependent differential form, is one of the basic equations that are studied in the sub-field of partial differential equations (PDEs), and has many applications in geometry, physics, scattering theory, integrable systems, etc.
In this work, after summarizing briefly (in Sect.2) some basic details of the Schrödinger equation, we focus on methods providing analytical solutions of this equation.Then (in Sect. 3), we present a numerical method based on artificial neural networks for solving the timeindependent Schrödinger equation.Next (in Sect.4), we apply our method in order to obtain numerical solutions for the wave functions describing the bound states of some exotic leptonic atoms and compare them with the corresponding analytic wave functions.Finally (in Sect.5), we summarize the main conclusions extracted from the present work and discuss future steps of our study.

Analytic solutions of the Schrödinger equation
For obtaining possible analytical solutions of the Schrödinger equation, the well known method of the separation of variables is utilized.For a system of two particles interacting via a potential V ( r), the time-dependent Schrödinger equation is formulated as [14,15], where = h/2π represents the well known reduced Plank constant.Assuming that the masses of the two particles are m 1 and m 2 (with m 1 ≥ m 2 ), then the reduced mass of the system, µ, entering Eq. ( 1) is defined as Notice that, for m 1 >> m 2 the resulting reduced mass is µ ≈ m 2 .
Then, if the potential V ( r) is spherically symmetric, i.e. if it holds V ( r) ≡ V (r), we apply again the separation of variables method to separate the radial wave function R(r) from the angular part of ψ( r) (here we take for granted that the angular dependent part is a spherical harmonic of the type Y ,m (θ, φ)), as ψ( r) = R(r)Y ,m (θ, φ).This procedure, results to a radial differential equation of the form We may, furthermore, simplify the latter equation by using the transformation R(r) = u(r)/r and obtain the reduced radial differential equation, In this work, we are interested in the solutions of the radial equation ( 3), in the special case for which V (r) is a Coulomb type potential given by (note here that the square of elementary charge is e 2 = 1.4399764MeV fm).In this specific case, the differential equation satisfied by the reduced radial wave function u(r) is written as The radial equation (3) admits analytic solutions of the form [15], where L 2 +1 n− −1 (x) denote the associated Laguerre polynomials which are defined as, The square root coefficient in Eq. ( 7) comes out of the normalization condition written as which must hold for the Radial wavefunction R(r) and the reduced radial wavefunction u(r).
As mentioned in the Introduction, the energy eigenvalues are determined simultaneously with the corresponding wave functions and are given by We note that only the quantum number n (primary quantum number) enters the definition of E. In Eq. ( 9), a 0 denotes the Bohr radius of the particle of mass m 2 orbiting around the attractive particle of mass m 1 placed at the position r = 0 (in the relative coordinate system) and is defined by where α is the well known fine structure constant In the latter definition 0 stands for the "vacuum electric permittivity"constant (c the speed of light).

Numerical solution of Schrödinger equation within N.N. techniques
Schrödinger's equation doesn't always posses analytic solutions, consequently a numerical approach must necessarily be employed.In this Section we present an advanced method derived within the context of neural networks (NN) which solves numerically the time independent Schrödinger equation ( 6) using the method developed by Kosmas and Lagaris [8].
Assuming that u(r) is a trial solution of Eq. ( 6), we proceed with the following steps.At first, we define a grid from r = 0 up to a point r = b, where the wavefunction is practically vanishing (see below Sect.4), and denote it by r i , for i = 1, 2, . . ., s.Then, because equation ( 6) must hold at every point r i of the grid, we write, equivalently, the expression IC-MSQUARE-2023 Journal of Physics: Conference Series 2701 (2024) 012133 The next step towards solving the latter equation and obtain the u(r), consists in (i) parametrizing u(r), and (ii) minimizing the left-hand side of equation (12).In order to avoid the trivial solution u(r) = 0 everywhere in the domain of r, we divide the above expression by the normalization |u(r)| 2 dr.To facilitate the minimization process, we use a parametrization of the form u(r) = re −kr N (r, u, v, w) where k ∈ R is an optimization constant (in our case k ≈ −2Eµc 2 / c), u j , v j , w j ∈ R d are the parameters which need to be determined through the minimization process.
The neural network N (r, u, v, w) entering Eq. ( 13) is defined as where f (s) denotes the activation functions employed.
In our present work, we utilized the following three activation functions [16]: (i) The sigmoid activation function defined as (ii) The ReLU activation function defined as (ii) The Softplus activation function defined as It is worth noting that, the function N (r, u, v, w) of Eq. ( 14) constitutes a 2-layer neural network.The reason we parameterized u in this way is because it satisfies the conditions u(0) = 0 and u(r) → 0 as r → ∞.By inserting u into the left-hand side of Eq. ( 12), we derive the error (or loss) function as Our task, in the last step, is to minimize the error function E f so as to approximate as close as possible the value zero.In practice, a value E f ≈ 10 −13 − 10 −10 is considered zero.We mention that, the effectiveness of the minimization becomes obvious from the comparison of the numerical function u(r) with the corresponding analytic solution (see below Sect.4).
In this work, we used two methods for the minimization procedure (in Python language): (i) the BFGS, and (ii) the trust-constr, both provided by the Python submodule scipy.optimize[17].Up to now, from these methods the performance of trust-constr proved to be better even though it requires longer training.
Regarding the determination of the energy E, this is calculated from the expression as long as the wavefunction u(r) is obtained.

Bound state wave functions for leptonic atoms
In this section, we use the numerical method derived in Sect. 3 in order to calculate the bound energy-spectra and the corresponding state wave functions in the purely leptonic (exotic) atom: Muonium (µ + , e − ).The results obtained for the radial wave functions R n, (r) for small quantum numbers n and are plotted in Fig. 1.
For the estimation of the level of performance we compare the numerical wave functions with the analytic ones given by Eq.( 7).Numerical results show clearly that our method is a very powerful numerical tool for solving the time independent Schrödinger equation.As can be seen from Fig. 1, the performance of our numerical method is reduced as n gets larger.Regarding the numerical solution of these  [9,10] utilizing different methods such as the Simulated Annealing [12], etc.
In Table 1, the corresponding values of the energies, obtained by using the eigenvalue formula (9) and the integral formula (??energint) are listed.We see that the agreement between the analytical and numerical results is very good.

Conclusions -Outlook
In this work, we derived a numerical method (in Python language) for solving the time independent Schrödinger equation.The algorithm is relying upon a two-layer artificial neural network.Then, we applied this method in order to solve this differential equation and determine the low-lying bound states (for n = 1, = 0 and for n = 2, = 0 and = 1) in purely leptonic atoms focusing on the Muonium (µ + ,e − ) leptonic atom.
From the calculations performed in this work, we concluded that for larger values of the principal quantum number n and > 0 a better training is required in order to achieve high confidence level of the agreement between the analytical and numerical wave functions.
In the next step, are going to extent the the Dirac equation since it is a relativistic description.

Table 1 .
Analytical and Numerical energies of the low-lying bound states of Muonium (µ + , e − ), purely