Confined hydrogen atom: endohedrals H@C36 and H@C60

In this work, for the lowest states with angular momentum, l=0,1,2 the energies and eigenfunctions of the endohedrals H@C36 and H@C60 are presented. The confining spherically-symmetric barrier was modeled by an inverted Gaussian function of depth ω 0, width σ and centered at r c , w(r)=−ω0exp[−(r−rc)2/σ2] . The spectra of the system as a function of the parameters ( ω0,σ,rc ) is calculated using three distinct numerical methods: (i) Lagrange-mesh method, (ii) fourth order finite differences and (iii) the finite element method. Concrete results with not less than 11 significant figures are displayed. Also, within the Lagrange-mesh approach the corresponding eigenfunctions and the expectation value of r for the first six states of s,p , and d symmetries, respectively, are presented as well. Our accurate energies are taken as initial data to train an artificial neural network that generates faster and efficient numerical interpolation. The present numerical results improve and extend those reported in the literature.


Introduction
In the majority of problems the time-independent Schrödinger equation, the central object in quantum mechanics does not admit exact solutions in terms of elementary or special functions. In such cases, it is necessary to solve this equation by means of numerical and approximate methods which inherently carry a certain degree of accuracy. For instance, some of the most common numerical methods are the following: the Numerov method, the spline-based method, the finite difference, the finite element, the Lagrange mesh method, and the variational one, to name a few. Interestingly, artificial intelligence algorithms have emerged as a promising tool for tackling eigenvalue problems of differential and integrodifferential operators [1,2] with important applications on the solution of Schrödinger equation [3][4][5][6][7] and open quantum systems [8][9][10]. In particular, artificial neural networks have shown to be accurate and efficient algorithms for calculating eigenvalues and eigenfunctions, due mainly to their high capabilities to approximate complex functions from data [11][12][13].
Classical numerical schemes are widely used in the study of both spatially confined and unconfined quantum systems. In particular, the investigation of spatially confined quantum systems has gained much interest because some of their physical properties change abruptly with the size of the confining barrier. Furthermore, many physical phenomena can be modeled by a confined quantum system for example: atoms and molecules subject to high external pressures, atoms and molecules in fullerenes, inside cavities such as zeolite molecular sieves or in solvent environments, the specific heat of a crystalline solid under high pressure, etc. A complete list of applications can be found in the articles and reviews on the subject [14][15][16][17][18][19][20][21][22][23][24][25][26].
Eighty-five years ago, Michels et al [27] calculated the variation of the polarizability of hydrogen under high pressure. To this end, they proposed a simple model of the confined hydrogen atom where, to a first approximation, the nucleus is anchored in the center of a spherical box of radius r 0 and impenetrable walls. It is assumed that this infinite potential is due to the presence of neighboring negative electric charges. The corresponding wave function vanishes at the surface of the sphere, i.e. it must obey Dirichlet boundary conditions. Since then, this model has been successfully applied in the study of confined many-electron atoms and molecules (see [28,29] and references therein). However, it takes into account the effects of repulsive forces only. A more realistic potential that embodies an attractive force considers a softer confinement in cavities of penetrable walls. The simplest penetrable confining potential of this type is that of a step potential [30,31].
Another relevant example of a penetrable potential is realized by the logistic potential used in the analysis of the entropic properties of hydrogen atom [32].
Likewise, in the study of atoms inside fullerenes, the corresponding potential has been modeled effectively by the attractive spherical shell [18], the δ-potential [33] and a Gaussian spherical shell [34,35]. Moreover, a spatially confined system inside an impenetrable spherical barrier, in addition to the inverted Gaussian potential [36][37][38][39], was employed by Lumb et al to describe the emission of two photons in the hydrogen atom [40].
It has been recognized as an advantage to study free systems, where bound states exist and the potential vanishes V → 0 at large distances, using the methods applied to the spatially confined systems (see [41][42][43]). In such a procedure the original system with potential V is additionally enclosed inside an impenetrable box of radius r 0 . Afterward, the wave functions as well as the eigenvalues are calculated in a sequence by increasing the value of this parameter r 0 until the desired accuracy is achieved. Mathematically, it has been demonstrated that this procedure converges to the exact results of the original system where the spatial confinement is absent (see [41][42][43]). Hence, the artificial spatial confinement gives highly accurate results with respect to the exact values, which can not be guaranteed by other methods [41]. Some examples to illustrate this point are the free hydrogen atom, an electron in the Yukawa potential, the Hulthen potential, and the singular oscillator (see [41][42][43]).
In this work, in the Born-Oppenheimer approximation, we calculate the energies and wave functions of the lowest states of a hydrogen atom confined inside a spherical penetrable barrier. This system has been used in [34,35] (and references therein) to model the hydrogen atom confined by C 36 and C 60 fullerenes, where the corresponding confining potential is described by an inverted Gaussian function. Here, a more accurate and systematic study is carried out. To solve the Schrödinger equation we employ three different numerical methods: the Lagrange mesh, finite difference, and finite element. They are complemented by the training and use of an artificial neural network. For the sake of self-consistency, a brief description of the aforementioned numerical methods is presented in the next section with the corresponding references.
The goal of the present study is threefold. Firstly, we aim to improve previous results reported in the literature. Secondly, we make a comparison of the numerical results obtained with different methods and discuss the advantages of each individual approach. Finally, we show that the complementary use of an artificial neural network leads to efficient numerical interpolation.

Methodology
The Schrödinger equation of hydrogen atom confined by a C 36 and C 60 endohedral cavities, in atomic units (ℏ = m e = e = 1), is of the form where ∆ is the three-dimensional Laplacian, and the Gaussian spherical barrier w(r) is given by here, ω 0 is the well depth, r c is the position of the center of the peak and σ is the width of the Gaussian, respectively. In concrete calculations, for the purposes of comparison with the existing results in the literature, from time to time the parameters r c and σ will be presented in Angstroms (Å). As is usual for any central potential, the angular momentum is conserved and the solutions of the Schrödinger equation (1) can be factorized in spherical coordinates (r, θ, ϕ). Explicitly, they read with Y l,m (θ, ϕ) being a spherical harmonic function. In the standard approach, it is convenient to introduce an auxiliary function u nl (r) defined by the relation R nl (r) = u nl (r)/r. In this case, from (1) and (3) it follows that u nl (r) must obey the isospectral radial problem  with the boundary condition u nl (r = 0) = 0. The spectral problem (4), is defined in the domain (configuration space) r ∈ [0, ∞). The effective potential appearing in (4) is given by l = 0, 1, 2, . . . , (n − 1) is the quantum number of angular momentum. Formally, (4) describes a one-dimensional particle of unit mass in the half positive line with an effective potential V eff (r). The shape of V eff is displayed in figure 1 for fixed l = 0, ω 0 = 1 a.u., σ = 0.5 a.u. and four different values of r c (the position of the center of the peak); there exists a critical value r crit c ≈ 1.05 a.u. such that for r c > r crit c the effective potential develops a second minima. In turn, figure 2 shows the behavior of V eff (5) as a function of the width σ of the Gaussian for fixed ω = 0.5 a.u. and r c = 3.54 Å.
In order to solve the radial Schrödinger equation (4), three different methods will be employed: (i) The Lagrange-mesh method (LMM), (ii) finite difference method (FDM), and (iii) finite element method (FEM). By combining our accurate results with an artificial neural network, we construct an efficient numerical interpolation for the energies.
at the N mesh points are given by The coefficients λ i are the weights associated with a Gauss quadraturê In terms of the N Lagrange functions f i (x) (7), the solution of the Schrödinger equation (4) is expressed as The trial function (9), together with the Gauss quadrature (8) and the Lagrange conditions (6) leads to the system of variational equations where ω(x i ) is the Gaussian potential (2) evaluated at the mesh points x i and T ij are the kinetic-energy matrix elements whose explicit expression is found in [46]. h is a scaling factor that allows to adjust the mesh to the system in consideration. By solving the system (10), not only the energies E are obtained but also the eigenvectors c i from which the approximation to the wave function (9) is obtained. Inside the LMM approach, the expectation value of the radial coordinate r is easily calculated. Given the approximation to the wave function (9) together with Gauss quadrature (8) and the Lagrange condition (6) leads to where c i are the eigenvectors resulting from solving (9) and x i are the mesh points.

Finite difference method
The FDM is a numerical approach easy to implement on a computer, and it is used to solve ordinary as well as partial differential equations in an approximate way. The method is based on the discretization of the Hamiltonian on a spatial grid. Consider the system defined in a domain [0, L], the Schrödinger equation is then solved on a uniform grid [47] defined by the set of discrete points which are the nodal points { r j } . The domain is split into N subintervals of equal length h: where In the present work, we use a fourth-order centered difference approximation to the second derivative, namelyũ Substituting this definition on the Schrödinger equation (4) we find that there is one equation for each j value. For j = 1 and j = N we impose the boundary condition ψ(r 1 ) = ψ(r N+1 = L) = 0. Hence, these systems of equations can be written as an eigenvalue problem H C = E C where the matrix H is a pentagonal matrix, whereas the vector C contains the values ofũ evaluated on the grid, i.e. C i =ũ(r i ). To solve the differential equation, we select the parameters of the potential and angular momentum l of the states, we propose a value of L much higher than r c , and solve the eigenvalue problem. We increase the value of L and solve the eigenvalue problem again, compare the energy values and repeat this process until the desired accuracy is achieved. We call L max (l) to the value of L at which the desired accuracy is achieved. Note that the boundary condition at ψ(r N+1 = L) = 0 is equivalent to the existence of an impenetrable barrier at that precise point, this fact establishes the relationship with the confined systems. In this case, the Schrödinger equation to be solved can be written as (cf (4)): where V eff (r) is the effective potential (5) and the confinement potential V c is defined as follows The above spectral problem (15) is defined in the reduced domain r ∈ [0, L]. The Schrödinger equation (15) was studied by Lumb et al [40] in connection with the two-photon emission in confined hydrogen atom. In the region r < L the Schrödinger equation becomes cf (4). The functionũ satisfies the Dirichlet boundary conditions: The original energy eigenvalues and eigenfunctions of the free (spatially unbounded) system (4) are recovered in the limit L → ∞,ũ

Finite element method
The FEM combines the discretization of the space in elements (intervals) and the use of polynomial interpolating functions on each element. This method is widely used in engineering, classical physics, and quantum mechanics problems, among others. The discretization is based on the reformulation of the relevant differential equation as an equivalent variational problem. The Garlekin methods are then employed in the corresponding minimization procedure [48]. In general, one can identify the following steps to solve the associated differential equation: (i) to present the problem in a variational formulation, (ii) a discretization of the domain using FEM, and finally, (iii) to find the solution of the discrete problem, which may consist of the solution of a system of simultaneous equations or an eigenvalue problem. In Quantum Mechanics the FEM was used from a few years ago [47,[49][50][51][52][53][54]. An excellent introduction to FEM in Quantum Mechanics can be found in the Ram-Mohan's book [52]. Only the key points of the method will be presented here.
The time independent Schrödinger equation for a particle of mass m subjected to a potential energy V(r) is given by: This equation can be obtained by minimizing the following action integral I: where ψ * and ψ are considered as two independent 'fields' . It is assumed that ψ is continuous up to its second derivative. By varying the action I with respect to ψ * we obtain the Schrödinger equation (20). For a spherical symmetric potential V = V(r), only the dynamics of the r coordinate is not trivial. The problem reduces to the one-dimensional case in the r-space, the variable r varies in the interval [0, L]. Hence, the action integral I (in atomic units) becomes: Now, the interval [0, L] is divided into small subintervals called elements. The action integral (22) can be decomposed as the sum of the action computed in each element, where n is the number of elements and I ( j) is the action integral evaluated on the jth element. Explicitly, the wave function ψ j defined in the jth element is expanded as a linear combination where c j , j = 1, 2, . . . , n, are unknown coefficients to be determined whereas N j (r) are interpolating polynomials. These are defined for the jth element and they are identically zero out of this element. The basic idea is to make the variation of the action integral with respect to the coefficients c * i , and by solving these equations to obtain the optimal energies and eigenfunctions.
In particular, Guimaraes and Prudente [53] developed an alternative version of FEM called the p-finite element method (pFEM) to study the confined hydrogen atom, and Nascimento et al [34] employed a pFEM version to study the electron structure of endohedrally confined atoms using an attractive Gaussian potential to model atoms inside fullerenes.
As L becomes very large the energies of a confined system approach to those of the confinement-free system. In practice, a large value of both L and the parameter n, and a polynomial degree for the N j (r) appearing in (24) are chosen and then, the generalized eigenvalue problem is solved. For a fixed L, by increasing either the value of n or the polynomial degree, or both, higher precision in the results can be achieved as we will explain in the next section.

Artificial neural networks
Artificial intelligence has emerged as a collection of computational techniques which seek to mimic the human brain in order to complete tasks for which standard algorithms lead to partially satisfactory results or are costly to implement [55][56][57][58]. Particularly, neural networks are artificial intelligence algorithms inspired by the workings of neurons in the human brain. These algorithms have demonstrated the capacity of pinpointing relevant specific pieces of information 'buried' in huge data sets and unveiling complex non-linear relationships between the inputs and target, which would be all but impossible to accomplish through a standard visual inspection [59][60][61].
In this work, we implement a neural network to estimate the eigenvalues of the radial Schrödinger equation (4) for different positions of the peak of the Gaussian, r c . This neural network consists of a two-layer feed-forward architecture with ten neurons in the hidden layer and a linear neuron in the output layer. Note that the input node of our neural network corresponds to the center of the inverted Gaussian potential, r c . Figure 3 shows the architecture of the regression neural network. In general, the output of each neuron before the activation function reads, where ω i are the synaptic weights, x i are the inputs, and N is the number of inputs. Importantly, all the neurons in the output layer contain linear activation functions whereas the neurons in the hidden layer have sigmoid functions given by Synaptic weights of the neural network are optimized with the Levenberg-Marquardt backpropagation method in a direction that minimizes the mean squared error [62,63]. This method approaches second-order training speed without having to compute the Hessian matrix. Because the performance function is given by a sum of squares then the Hessian matrix can be approximated as H = J T J, where J is the Jacobian matrix. Using this approximation, the synaptic weights can be updated by the following expression, where k denotes the kth iteration, µ is the learning rate and e is the vector of networks errors. To train the neural network, we use a subset of the dataset that contains eigenvalues of the radial Schrödinger equation (4) for different positions r c of the peak of the Gaussian. These eigenvalues are calculated by the LMM. After the training stage, the neural network is able to predict the eigenvalues of the whole dataset with a coincidence at six significant digits at least.

Results and discussion
In order to describe the confinement of the hydrogen atom inside C 36 and C 60 fullerene cages, an attractive spherical Gaussian potential [34] is considered (1). Energies, eigenfunctions and expectation values of r for the first six states of s, p and d symmetries are accurately calculated. In order to be able to compare with previous results, the values of the parameters of the Gaussian potential (1) that are investigated in detail are: r c = 2.5 Å (for H@C 36 ) and r c = 3.54 Å (for H@C 60 ), ω 0 = 0.5 a.u. and σ = 0.26, 0.57 and 1.59 Å. Previous work has shown [64][65][66] that the Lagrange-mesh (LMM) method provides very accurate values for the energy values. The corresponding results obtained using the LMM are shown in tables 1-3 with 12 significant figures; overall within this precision, as a result of explicit computations it turns out that the LMM is the most practical, accurate and easier to implement an approach for the present case (see discussion below). We emphasize that, in principle, the three distinct numerical approaches (i) LMM, (ii) fourth order finite differences, and (iii) the FEM, can lead to extremely highly accurate results by choosing more refined values of the parameters of the methods with the corresponding growth of the computation time.
Before discussing the concrete results let us briefly mention some details about the LMM. In the system (10), there are two free parameters: the size of the mesh N and the scaling factor h. The optimal values of these two parameters (N, h) depend on the considered state of the system as well as on the value of the parameters occurring in the Gaussian potential. In general, the results presented in tables 1-3 for states s, p and d, respectively, are obtained with a mesh of at least N = 250 points and h in some interval between h ∈ [0.02 − 1.0]. For states with small quantum numbers, even meshes of size n = 120 are adequate. For example, if we consider the state 1s with w 0 = 0.5, r c = 2.5 Å and σ = 2.26 Å the energy is obtained by using (N, h) = (120, 0.04 − 0.12) while for the state 8d, for the same values of the parameters of the potential, (N, h) = (250, 0.26 − 0.30). Table 1. Energy E and expectation value ⟨r⟩ of the hydrogen atom of the 1s . . . 6s states in presence of a inverted Gaussian potential (2) centered at rc = 2.50 and 3.54 Å with w0 = 0.5 a.u. as a function of σ. The free case corresponds to σ = 0. For rc = 3.54 Å comparison is done with results presented in [35]. The values for E and ⟨r⟩ are displayed with 12 significant figures, and they correspond to the results obtained with the Lagrange Mesh Method. The coincidence with the FEM results is in 7-8 digits whilst for the FDM, the agreement with the results of the LMM is in all figures. The convergence of the LMM is determined by the stability of the results with respect to an increase in the size of the base N and the variation of h. Table 1 presents the results for the 1s, . . . , 6s states. For r c = 3.54 Å a comparison with [35] is possible for the levels 1s, 2s, 3s, 4s with σ = 1.59, 0.57 and 0.26 Å where it can be seen a complete agreement in 7 decimal digits (for 1s and 2s states) and 8 decimal digits (for 3s and 4s). For completeness, the case σ = 0 (the free hydrogen atom) is also presented. For all these s-states, the energy decreases by increasing the value of σ as can be seen in figure 4(c). On the other hand, the presence of the Gaussian potential has an important effect on the expectation value of r as indicated in the seventh column of table 1 (see also figure 4(d)): at fixed r c = 3.54 Å (i) as a function of σ ∈ [0.00, 1.59] Å, ⟨r⟩ increases for the ground state, whilst (ii) for the states 4s, . . . , 6s, ⟨r⟩ decreases. Columns 3 and 4 of table 1 (see also figures 4(a) and (b)) present the results of the energy and the expectation value ⟨r⟩ when the center of the Gaussian potential is r c = 2.50 Å. As well as for r c = 3.54 Å, the energy decreases as a function of σ. The behavior of the expectation value ⟨r⟩ is depicted in figure 4(b). This effect on ⟨r⟩ reflects how the electronic charge is attracted to the Gaussian part of the potential.
States with l = 1 (2p, . . . , 7p) are presented in table 2 for r c = 2.50 and 3.54 Å. In both cases, with respect to the free configuration (σ = 0), the system gets more bound as the value of σ increases (see figures 5(a) and (c)). For r c = 3.54 Å a comparison with [35] is also possible: we see a complete agreement in 7 decimal digits for the 2p state and 8 decimal digits for the 3p and 4p states for the three values of σ = 0.26, 0.57 and 1.59 Å. The expectation value ⟨r⟩ is shown in figures 5(b) and (d).
Results of the energy and the expectation value ⟨r⟩ for l = 2 (3d, . . . , 8d) are displayed in table 3 for r c = 2.50 and 3.54 Å. For both values of r c , by increasing σ the energy becomes more negative compared to the unconfined system (σ = 0). The expectation value ⟨r⟩ exhibits a decreasing behavior with the increasing of σ for all states except the lowest 3d state, which decreases and eventually increases.

Comparison of LMM with the finite difference and FEMs
The FDM was implemented in the second and fourth order in Matlab. Calculations using second-order finite differences require a very large number n of nodal points, which implies to diagonalize very large matrices, and therefore a significant computational time is involved. For this reason, we decided to use fourth-order Table 2. Energy E and expectation value ⟨r⟩ of the hydrogen atom of the 2p, . . . , 7p states in presence of a inverted Gaussian potential (2) centered at rc = 2.50 and 3.54 Å with w0 = 0.5 a.u. as a function of σ. The free case corresponds to σ = 0. For rc = 3.54 Å comparison is done with results presented in [35].     provides results such that the rounding error does not affect the first 11 decimal digits. Another quantity that can affect the accuracy of the eigenvalues is the quantum state to be calculated. For example, at fixed values r c = 3.54 Å, ω 0 = 0.5 a.u. and σ = 0.26 Å, in order to obtain eigenvalues with 8 decimal places of accuracy it is necessary that L max = 150 a.u. and a grid with n = 15 000 for states with l = 0, whereas for states with l = 1, L max = 170 a.u. and n = 17 000, and for l = 2 the optimal parameters are L max = 200 a.u. and n = 20 000. This shows that to calculate excited states, while maintaining the same accuracy of the previous energy levels, it is necessary to increase the value of L max as well as n. Additionally, by increasing the value of L max and n the agreement with the results of the LMM is in all digits. When the FEM is applied, the agreement with the results of the LMM is in 8 decimal digits. It should be mentioned that the MATHEMATICA 12.3 software package allows us to calculate the solutions of the eigenvalue problem (4) as well. It can be easily done using the NDEigensystem command (based on the FEM) which, in general, provides the smallest eigenvalues and eigenfunctions of the involved linear differential operator on a certain finite region, L max ∼ 150. For the lowest states, in table 4 we display the relative difference |EM−ELM| |ELM| between the energy E M obtained with MATHEMATICA 12.3 and the corresponding value E LM computed in the Lagrange-Mesh approach. Specifically, at fixed values of s-states (l = 0) with parameters r c = 3.54, ω 0 = 0.5, and σ = 1.59, by using a maximum cell measure of 0.001 the optimal value of L max increases monotonically from L max = 30 at n = 1 3 (ground state) up to L max = 120 for n = 6. The calculations were run in MATHEMATICA 12.3 on a personal laptop.

Numerical interpolation via an artificial neural network
Additionally, we trained a regression neural network to estimate the energies of the hydrogen atom for the 1s state as a function of the center r c of the inverted Gaussian potential (2). It is worth mentioning that the neural networks were trained and tested using the Neural Network Toolbox from Matlab 2022a, which runs on a computer with an Intel Core i7-2600 CPU (@3.40 GHz) and 16 GB of RAM. The training and testing data are generated by the LMM. As is standard in artificial neural networks, we devoted 70% of the dataset to training, 10% to validation, and 20% to testing. After training, the neural network can predict the Lagrange-mesh results with not less than six significant digits which is the relevant domain in the present non-relativistic framework. Despite we show solely the results for the Lagrange method in table 5, we also performed the training of neural networks using the FEM and FDM data, respectively. In both cases, the obtained results present the same six significant digits as those obtained with the Lagrange mesh method. A plausible explanation is that all methods provide accurate energies with not less than 8 significant figures, which always is over the accuracy that neural networks can solve for the dataset. Remarkably, our algorithm takes 40 µs to calculate the energy for a given value of r c . For ω 0 = 0.5 a.u. and σ = 0.26 Å, table 5 displays the results obtained for the energy as a function of r c expressed as r c = λ r 0 where r 0 = 3.54 Å and λ is a factor specified in the first column. As can be appreciated in table 5, there exists a critical value r c ≈ 1.011 a.u. for which the energy of the system reaches its minimum value. In general, for excited states, the energy presents maxima and minima as a function of r c .
In the present computations, the depth of the well was chosen as ω 0 = 0.5 a.u.. For larger values of this parameter, the confined hydrogen atom energy becomes more negative in complete agreement with the results reported in [34]. When possible, in each of the tables, the obtained results in the present work are compared with those presented by Lin and Ho [35]. As can be seen, the methods used in our study lead to an improvement in the energy values reported previously.
In the sections where the results of each specific methodology are presented, the calculation conditions are indicated. However, the reader may contact the authors for further information. Table 5. Energies of the confined hydrogen atom for the 1s state as a function of the center of the inverted potential (2), rc. Here, ELM is the energy calculated with the LMM and EAI is the energy with the trained neural network using as an input the Lagrange-Mesh values. These results were obtained using the following parameters of the potential: ω0 = 0.5 a.u., σ = 0.4913287924027 a.u. and rc = λ r0 with r0 = 6.6896304811752 a.u.

Conclusions
In the Born-Oppenheimer approximation, for the lowest states with angular momentum l = 0, 1, 2 the energies and eigenfunctions of the hydrogen atom inside a C 36 and C 60 fullerenes are presented. The confining barrier was modeled by an inverted Gaussian function w(r) = − ω 0 exp[−(r − r c ) 2 /σ 2 ]. The approximate solutions of the corresponding Schrödinger equation were determined by three different numerical methods: (i) the LMM, (ii) the (fourth-order) finite difference and (iii) the FEM. As a complementary tool, we use an artificial neural network to interpolate/extrapolate the results. In general, the confined hydrogen atom becomes more bound with the increase of the width σ and the depth ω 0 of the Gaussian potential, for the two considered cases of r c . Using the LMM accurate energies with not less than 11 significant figures were obtained. The optimal values of the size of the mesh N and the scaling factor h depend on the state being studied as well as on the parameters of the confining potential w(r). Since this method is not completely based on a variational principle, we must point out that the computed energies are not necessarily greater than or equal to the (unknown) exact ones. However, in all known cases where the Lagrange mesh results are stable with respect to the variation of N and h, like in the present study, it turns out that it converges rapidly and generates simple highly accurate solutions.
The FDM is a suitable method for solving the radial Schrödinger equation (14). In the present work, we used a fourth-degree approximation for the kinetic energy operator (18). The accuracy of the eigenvalues depends strongly on the corresponding step size h (equation (19)). In particular, using a fourth-degree approximation the error in the kinetic energy operator (equation (18)) is of order h 4 . The calculations were performed in a code written with double precision in MATLAB. It should also be noted that this method, like the LMM, is not based on the variational principle.
In the case of the FEM, precise energies (always from above the exact ones) can be calculated by simultaneously increasing L max , the number of elements in which the interval [0, L max ] is divided and the degree of the interpolating polynomials. Moreover, this method allows us to deal with problems in higher spatial dimensions with regular and irregular boundaries, which is not so easy to implement in the Lagrange-mesh and FDMs.
Finally, the artificial neural network is computationally a faster efficient tool to compute the spectra. Nevertheless, for the training stage, it requires to know in advance accurate results in several points on the space of parameters. In combination with the Lagrange-mesh or the FDM it significantly reduces the overall computational time, although with less accuracy. By means of the methods used in the present study, energies were obtained with higher accuracy than those reported in the literature. In future works, we plan to establish in a more systematic manner the implementation of artificial neural networks in the study of other physically relevant atomic and molecular systems.

Data availability statement
No new data were created or analysed in this study.