Fisher Information-Based Optimization of Mapped Fourier Grid Methods

Abstract

The mapped Fourier grid method (mapped-FGM) is a simple and efficient discrete variable representation (DVR) numerical technique for solving atomic radial Schrödinger differential equations. It is set up on equidistant grid points, and the mapping, a suitable coordinate transformation to the radial variable, deals with the potential energy peculiarities that are incompatible with constant step grids. For a given constrained number of grid points, classical phase space and semiclassical arguments help in selecting the mapping function and the maximum radial extension, while the energy does not generally exhibit a variational extremization trend. In this work, optimal computational parameters and mapping quality are alternatively assessed using the extremization of (coordinate and momentum) Fisher information. A benchmark system (hydrogen atom) is employed, where energy eigenvalues and Fisher information are traced in a standard convergence procedure. High-precision energy eigenvalues exhibit a correlation with the extrema of Fisher information measures. Highly efficient mapping schemes (sometimes classically counterintuitive) also stand out with these measures. Same trends are evidenced in the solution of Dalgarno–Lewis equations, i.e., inhomogeneous counterparts of the radial Schrödinger equation occurring in perturbation theory. A detailed analysis of the results, implications on more complex single valence electron Hamiltonians, and future extensions are also included.

1. Introduction

The quantum mechanical treatment of practically any problem in atomic and molecular physics requires the solution of the time-independent Schrödinger differential equation to obtain energy eigenvalues and eigenfunctions. Furthermore, there are numerous situations where applying the Rayleigh–Schrödinger perturbation theory involves solutions for the inhomogeneous counterparts of the Schrödinger equation, with these counterparts characterized by the presence of an appropriate source term. These differential equations have analytical solutions only for a few specific simple cases; their majority, however, must be solved numerically. Evidently, the implicated numerical methods need to be validated for their efficiency and accuracy, and this is typically achieved by comparing numerical results to analytical solutions of benchmark problems. Still, a successful comparison does not automatically guarantee an equally efficient performance when exact solutions are not available. Hence, appropriate and robust convergence criteria are required that probe the efficiency and accuracy of a given numerical approach under any circumstances, particularly when variational methods are difficult to apply or not applicable at all.

Motivated by the above arguments, the purpose of the present article is to identify such convergence criteria. Specifically, our case study concerns atomic Hamiltonians, referring to a single valence electron subjected to a model local radial potential energy that exhibits a long-range Coulomb tail. The hydrogenic potential energy is naturally our benchmark system, but the approach can also be extended to more complicated cases as well.

There is a vast amount of earlier work involving the numerical solution of the Schrödinger equation in problems dealing with central potentials [1,2] or being part of computational methods like the convergent close-coupling [3] and the wave-packet convergent close-coupling [4] ones, the time-dependent close-coupling method [5], or the solution of the time-dependent Schrödinger equation with constant/variable-step spatial grids [6]. In the present work, for the direct numerical integration of the Schrödinger equation (with or without a source term), our method of choice is the grid-based discrete variable representation (DVR) where the basis functions are eigenfunctions of the coordinate operators [7,8,9,10]. Nowadays, this method and its variants are continuously refined and have proven to be quite efficient, either as stand-alone versions or in conjunction with finite element methods [11,12,13]. In any DVR, local operators depending on a given coordinate (such as the radial potentials) are diagonal in the coordinate basis. This property, however, emerges from the approximate treatment of the integrals defining the relevant matrix elements, and this simplification leads to the loss of the variational character of the calculated eigenvalues and eigenvectors [7]. It is only in the limit of an infinite dimension basis that convergence to the exact energy levels is expected. In fact, this lack of variational character when a finite dimension basis is employed accentuates the need to devise efficient convergence criteria, thus motivating us to use DVR as a Schrödinger equation solver. Along the same spirit, in the present work, we use a particular DVR variance, namely, the Fourier grid method (FGM), characterized by its exceptional simplicity. FGM’s discretization constraint (constant step size) is the only limiting factor in employing the method in the accurate solution of the Schrödinger equation for the Coulomb radial potential (hydrogen atom), as well as for more complicated ones exhibiting Coulomb tails at large distances. As duly noted by several authors, this deficiency can be successfully tackled by an appropriate coordinate transformation (mapping scheme) [12,13,14,15,16,17,18,19,20,21,22,23,24], leading to the so-called mapped-FGM. The choice of the most suitable mapping function is most frequently based on semiclassical and (classical) phase-space arguments [14]. These predictions, however, are occasionally unsuccessful [12,17], and this will also be demonstrated in the present work.

The preceding discussion suggests that the development of new methodologies capable of evaluating the relative convergence offered by different mapping functions is in order. Clearly, it is important for the convergence to be attested without any prior reference to phase-space arguments but based solely on the estimation of the relative computational error. Within this framework, the ever-evolving information theory of quantum mechanical systems appears as one of the most suitable methodological tools. A central point within information theory is the definition of appropriate entropic quantities that measure the information theoretic content of atomic and molecular charge densities. Among these quantities, most frequently discussed are Shannon entropy and Fisher information measures that were used for the description of several physical phenomena [25,26,27] and for constructing a number of important quantum uncertainty relations [28,29,30,31]. In particular, Fisher information is the central variable of the principle of extreme physical information and leads to the derivation of the Schrödinger wave equation itself [25,32]. Among the other works, the behavior of both Shannon entropy and Fisher information was examined with respect to the atomic number [33] and the degree of confinement (see, for example, [34,35,36,37,38,39]). More important in the present work is the fact that the entropic character of these quantities, and mainly their extremization, may allow them to act as variational evaluation parameters. Shannon entropy is a global measure of particle density spreading, while Fisher information is sensitive to the density’s local rearrangements [40]. Shannon and related entropies have been used in the past as a consistent measure of the quality of approximate wavefunctions [41,42,43,44], but it is the local character of Fisher information that renders it a more useful tool in assessing the precision offered by a given mapping scheme. In fact, this was demonstrated in related earlier studies, where its rate of change as a function of perturbative order was found to perform well as a convergence indicator in semiclassical perturbative quantum defect expansions [45]. Here, we extend these notions and explore their applicability to the variational optimization of diagonalization procedures, particularly the mapped-FGM-based ones. The success of this endeavor will allow for complementing existing error evaluation methodologies with new tools having their origin in information theory.

The rest of the paper is organized as follows: In the Section 2, we describe the main methodological aspects of FGM, along with the semiclassical and phase-space arguments that usually serve as selection criteria for the mapping function in mapped-FGM. In the Section 3, we present the expressions employed in this work for the coordinate and momentum Fisher information measure of free single-valence-electron atoms. Additionally, for the coordinate Fisher information, we derive the expression appropriate for s-state atoms embedded in a static electric field. In the Section 4, we present a series of Fisher information measure-based convergence tests for atomic energy levels and static dipole polarizabilities computed by the mapped-FGM. These tests demonstrate the variational character of the Fisher information measure and its ability to quantify the performance of a given mapping function with respect to the number of grid points and the radial extension of the calculation. Finally, in the Section 5, we conclude and discuss possible continuations of the present work. Atomic units (a.u.) are used throughout the paper.

2. Mapped Fourier Grid Method

In the present work, we are interested in the solution of radial equations of the form,

$$\left[-{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^2}{d{r}^2}}+U_{eff}\left(r\right)-E\right]P\left(r\right)=w\left(r\right).$$

(1)

In Equation (1), the effective potential energy is given by

$$U_{eff}\left(r\right)=U\left(r\right)+{\displaystyle \frac{\ell \left(\ell +1\right)}{2r^2}},$$

(2)

where ℓ is the orbital angular momentum quantum number. This study investigates atomic potentials with a long-range Coulomb tail and focuses on the free (unconfined) hydrogen atom, where U(r) = −Z/r and the nuclear charge is Z = 1.

The homogeneous form of Equation (1) (w(r) = 0) classifies as the Schrödinger equation (SE), while the inhomogeneous form of Dalgarno–Lewis (DL) equations is encountered in various forms in perturbation theory applications [21,46]. In the latter case, w(r) is an appropriate driving source term.

2.1. Basic Elements of a DVR and the Fourier Grid Method

Details on DVR may be found in numerous publications (see, for example, Refs. [7,8,9,10]), so we provide only a brief exposition of the key concepts here. Any DVR is a discretized coordinate representation. For the radial coordinate basis, |r_i> initially adopted in this work, scalar potential operators U are approximated by diagonal matrices <r_i|U|r_j> ≈ U(r_i)δ_ij, while the kinetic energy operator, $\widehat{T}=-\frac{1}{2}\left(d^2/d{r}^2\right)$, is approximated by a truncated N × N discrete matrix,

$$T_{ij}=\langle r_i|\stackrel{⌢}{T}|r_j\rangle i,j=1,2,3,\dots ,N.$$

(3)

By expanding over an arbitrary complete basis set, Equation (3) writes,

$$T_{ij}={\displaystyle \sum _{m,n}\langle r_i|m\rangle \langle m\left|\widehat{T}\right|n\rangle }\langle n|r_j\rangle$$

(4)

where |m> and |n> are state vectors of the chosen set, the latter dictating the nature of the radial grid (linear or nonlinear). In FGM, a Fourier basis is selected that is associated with an even number N of plane waves in momentum space, hence corresponding with N discrete values of momentum. Then, T_ij writes as [7]:

$$T_{ii}={\displaystyle \frac{{\pi }^2}{\mathrm{\Delta }{r}^2}}{\displaystyle \frac{N^2+2}{6N^2}}$$

(5a)

$$T_{ij}={\left(-1\right)}^{i-j}{\displaystyle \frac{{\pi }^2}{\mathrm{\Delta }{r}^2}}{\displaystyle \frac{1}{N^2{\sin }^2\left[\left(i-j\right)\pi /N\right]}}, i\ne j$$

(5b)

where Δr is the radial step size between the equally spaced discrete data points (linear, radial grid). Equations (5a) and (5b) can be further simplified in the limit of an infinite number of grid points N→∞, yielding,

$$T_{ii}={\displaystyle \frac{{\pi }^2}{6\mathrm{\Delta }{r}^2}}$$

(6a)

$$T_{ij}={\left(-1\right)}^{i-j}{\displaystyle \frac{1}{\mathrm{\Delta }{r}^2}}{\displaystyle \frac{1}{{\left(i-j\right)}^2}}, i\ne j$$

(6b)

It was shown [10] that the above representations correspond to a particle-in-a-box coordinate basis with “box” size R = N∙Δr. It was also pointed out [10] that Equations (6a) and (6b) are strictly valid only in the (−∞,+∞) interval, while the representation derived for the radial interval, [0,+∞) in the N→∞ limit writes,

$$T_{ii}={\displaystyle \frac{{\pi }^2}{6\mathrm{\Delta }{r}^2}}\left[1-{\displaystyle \frac{3}{2{\pi }^2{i}^2}}\right]$$

(7a)

$$T_{ij}={\displaystyle \frac{1}{\mathrm{\Delta }{r}^2}}\left[{\displaystyle \frac{{\left(-1\right)}^{i-j}}{{\left(i-j\right)}^2}}-{\displaystyle \frac{{\left(-1\right)}^{i+j}}{{\left(i+j\right)}^2}}\right], i\ne j.$$

(7b)

Equation sets (6a), (6b), (7a) and (7b) can be thought of as “universal”, as they are independent of the basis used in their derivation [10]. Therefore, they may be used as a “black box”. It is worth mentioning that in the present work, Equations (6a) and (6b) performed slightly better than Equations (7a) and (7b) in computing bound state energy levels, dipole polarizabilities and Fisher information measures and were thus selected as the form of choice for the rest of this work. Irrespective of the version, however, one should keep in mind that the number of grid points, N, is associated with a momentum cut-off, p_max, while box size, R, is associated with an energy cut-off, E_max. When solving the SE, these two cut-offs force an upper limit to the number of energy eigenvalues that can be computed with the desired precision [12,14,16,17,47].

2.2. Mapping Considerations and Procedure

Although the Fisher information-based optimization proposed in the present work is expected to be operational for any DVR variant, the particular present choice of FGM over other variants (see, for example, [48]) is guided by the simplicity implied by the constant step size grid. However, equally spaced grids are adequate solely for short-range potentials; long-range ones would require a computationally expensive high value of N. Fattal et al. [14] thoroughly discussed the inefficiency of a constant step-size FGM even for the simple hydrogenic potential. The origin of the inefficiency is twofold: First, the Coulomb potential singularity for r→0 and, second, its long-range nature (momentum space singularity). The difficulties are more important for s-states (ℓ = 0) that require a dense grid near the origin, while a sparse grid is adequate at large r. Clearly, the constant radial step size Δr must be abandoned in favor of a variable one. Variability can be enforced either by a suitable non-FGM basis set in Equation (4) or by an appropriate change of variable (mapping). The latter option is advantageous because it preserves the pros of FGM as the equidistant grid locations are reserved for the new variable.

The above considerations led Fattal et al. [14] and Kokoouline et al. [16,19] to propose a WKB-based, optimal, potential-specific transformation following the variation of the de Broglie wavelength of a particle subjected to an effective potential U_eff(r). Specifically, they proposed a transformation of the form,

$$s\left(r\right)={\displaystyle \frac{1}{\pi }}{\displaystyle \underset{r_{\mathrm{o}}}{\overset{r}{\int }}{\left[2\left(U_{eff,\max }-U_{eff}\left(r^{\prime }\right)\right)\right]}^{1/2}}d{r}^{\prime }$$

(8)

where U_eff,max controls the value of the inner radius r_o that should be slightly smaller than the inner turning point corresponding to the maximum energy of interest. Equation (8) leads to analytic s(r) functions for either the non-relativistic [14] or the relativistic [15] Coulomb problem. For more complicated atomic single-particle potentials, the calculation of Equation (8) should, in principle, be carried out numerically. This complication, however, may be avoided by following the proposal of Kokoouline et al. [16,19], where a simple analytic reference potential U_eff,env is chosen that is similar in shape and envelopes U_eff (i.e., it always holds U_eff,env ≤ U_eff and the two curves do not cross each other). This leads to an analytically known function s(r) and a grid that is denser everywhere than the true optimal one. In an earlier work [21] dealing with long-range coulombic potentials, the ℓ = 0 hydrogenic potential was employed as the reference one. Then, r_o = 0 and s(r) = br^1/2, where the constant b depends on the nuclear charge Z. Since in atoms heavier than hydrogen, the nuclear charge becomes a function of r, a new dimensionless variable was subsequently introduced that minimizes this r-dependence, namely,

$$x\equiv {\displaystyle \frac{s(r)}{s_{\max }}}$$

(9)

with s_max = s(r_max) and r_max is the maximum radius employed in the calculation (of the order of 10¹–10³ atomic length units for atomic scale calculations). The transformation of Equation (9) implies a mapping scheme of the form,

$$r\left(x\right)=r_{\max }h\left(x\right)$$

(10)

with the mapping function

$$h\left(x\right)=x^2$$

(11)

and the variable x restricted to the interval [0,1]. In practice, to avoid the singularities of U_eff(r), Equation (10) is reformulated with a spatial shift to r(x) = (r_max − r_init)h(x) + r_init, with r_init = 10⁻⁷ a.u. Hence, the radial “box” size is R = r_max − r_init ≈ r_max, which corresponds to a cut-off energy E_max = −Z/r_max. Instead, the mapped-FGM box size is now always R_x = 1, and the constant step-size is given by,

$$\mathrm{\Delta }x={\displaystyle \frac{1}{N-1}}$$

(12a)

that is, the x-grid is written as,

$$x_i=\left(i-1\right)\mathrm{\Delta }x, i = 1,2,3,\dots ,N.$$

(12b)

It has been shown [21,49] that the scheme of Equations (10) and (11) deals quite efficiently with all ℓ ≥ 0, whether for pure Coulomb potentials or more complicated ones with a long-range Coulomb tail. To understand its efficiency, let us follow the classical phase-space reasoning of Fattal et al. [14]. Then, for ℓ = 0 and Z = 1, we have,

$${\displaystyle \frac{1}{2}}{p}_r^2-{\displaystyle \frac{1}{r}}=E_{\max }$$

(13)

and, therefore,

$$p_r=\pm {\left[2\left(-{\displaystyle \frac{1}{r_{\max }}}+{\displaystyle \frac{1}{r}}\right)\right]}^{1/2}$$

(14)

from which it follows that p_r,max→∞ as r→0. Consequently, there is no reasonable rectangular FGM area, A_box = 2r_maxp_r,max, able to enclose the full phase-space area A of Equation (14),

$$A=2{\displaystyle \underset{0}{\overset{r_{\max }}{\int }}{p}_{r}dr}=\pi {\left[2r_{\max }\right]}^{1/2},$$

(15)

even though A itself is finite. In other words, since the value of p_r,max involved in a calculation is controlled by N, the full enclosure of A within A_box is ensured only in the impractical limit N→∞. Assuming now a general mapping scheme h(x), it can be easily shown that the transformed momentum is written as,

$$p_x=\pm {\left[2r_{\max }\right]}^{1/2}\left|{\displaystyle \frac{dh\left(x\right)}{dx}}\right|\cdot {\left[{\displaystyle \frac{1-h\left(x\right)}{h\left(x\right)}}\right]}^{1/2}.$$

(16)

A meaningful transformation h(x) should reshape the phase space in order to avoid singularities and minimize the wasted phase space of the original radial coordinate representation. Therefore, the maximum value p_x,max of p_x should be finite and, according to the same classical phase-space arguments, as small as possible. The rectangular mapped-FGM area is now A_box = 2R_xp_x,max = 2p_x,max, and the mapping quality may be estimated by inspection of the ratio Λ≡A/A_box, which should ideally approach unity. The transformation of Equation (11) offers p_x,max = p_x(x = 0) = 2[r_max]^1/2 and a quite large ratio Λ = π/4 (≈78.5% phase-space coverage).

The above semiclassical and classical considerations provide a recipe for the selection of h(x) through Equations (8) and (16). Their efficiency criterion is the maximization of Λ. Transformations that can, in practice, be proven more efficient than the one defined by Equation (11) but that do not comply with Equation (8) and/or lead to lower values of Λ, need to be (and have been) selected through either other non-WKB methodologies [18,23,24] or by making an educated guess [12,17]. Furthermore, the above arguments provide no practical hint on the optimal (N,r_max) sets. It is simply expected that, for a given N, large r_max values would result in quite sparse grids and loss of accuracy, while values that are too small may lead to the system’s confinement. As for Ν itself, it is, of course, anticipated that more accurate results are to be expected as the N→∞ limit is approached and the variational character of the method is restored. Since very large N values are impractical and time consuming, optimal sets of N, r_max, as well as of other parameters that could be incorporated in h(x), are most frequently found by examining the mapped-FGM performance to a known system and for a given targeted quantity (energy level, polarizability, etc.). Here, the Fisher information measure is supplementarily investigated for its potential as a convergence criterion.

2.3. Transformation of the Differential Equations

The mapping function of Equation (11) will be our main working transformation. To preserve the generality of the transformation in the subsequent calculations, h(x) is left unspecified in this sub-section. Furthermore, although other intervals have appeared in the literature [14,23,24], we shall continue to restrict variable x within the [0,1] interval and the step-size to be given by Equation (12a). Hence, the only restriction imposed on h(x) is that h(0) = 0 and h(1) = 1.

Using Equation (10), the Jacobian of the transformation from r to x is,

$$J\equiv {\displaystyle \frac{dr}{dx}}=r_{\max }{\displaystyle \frac{dh\left(x\right)}{dx}}=r_{\max }j\left(x\right)$$

(17)

where we have defined the function j(x) > 0, which dictates the variation of the radial step Δr(x). It turns out [16,17,19,20,23,24] that by writing the radial wavefunction as,

$$P\left(r\left(x\right)\right)=c\cdot j^{1/2}\left(x\right)\cdot y\left(x\right)$$

(18)

(with c a constant) one can eliminate the first derivative (d/dx) from the transformed differential equation for the auxiliary function y(x), which acquires the form

$$r_{\max }^{-2}\left[-{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^2}{d{x}^2}}+{\displaystyle \frac{3}{8}}{\left({\displaystyle \frac{j^{\prime }}{j}}\right)}^2-{\displaystyle \frac{j^{″}}{4j}}+r_{\max }^2{j}^2\left(U_{eff}\left(r\left(x\right)\right)-E\right)\right]y=j^{3/2}w\left(r\left(x\right)\right)$$

(19)

The primes in Equation (19) denote differentiation with respect to x. Elimination of the first derivative is desirable to keep the transformed kinetic energy operator symmetric. Indeed, exploiting the aforementioned universality of FGM kinetic energy matrices, the matrix of the new “kinetic energy” operator ${\widehat{T}}_x=-\frac{1}{2}\left(d^2/d{x}^2\right)$ is simply given by any set of Equations (5)–(7) with the substitution Δr→Δx. Then, the solution of SE (w = 0) reduces to the generalized eigenvalue problem,

$$r_{\max }^{-2}Ay=EBy$$

(20a)

where the elements of the symmetric matrix A and the diagonal matrix B are given by

$$A_{ij}=T_{x,ij}+\left({\displaystyle \frac{3}{8}}{\left({\displaystyle \frac{j^{\prime }\left(x_i\right)}{j\left(x_i\right)}}\right)}^2-{\displaystyle \frac{j^{″}\left(x_i\right)}{4j\left(x_i\right)}}+r_{\max }^2{j}^2\left(x_i\right)U_{eff}\left(x_i\right)\right){\delta }_{ij}$$

(20b)

and

$$B_{ij}=j^2\left(x_i\right){\delta }_{ij}$$

(20c)

By defining another auxiliary function φ(x), such that y(x) = j^−1/2φ(x), Equation (20a) is transformed back to an ordinary eigenvalue problem [23,24],

$$r_{\max }^{-2}{B}^{-1/2}A{B}^{-1/2}\phi =E\phi$$

(20d)

where $B_{ij}^{-1/2}={\delta }_{ij}/j\left(x_i\right)$, and the matrix B^−1/2AB^−1/2 remains symmetric.

On the other hand, when dealing with DL equations (w ≠ 0), there is no need for such a y→φ transformation. Nevertheless, we adopt it here to unify our treatment. The DL equation reduces to the linear system

$$r_{\max }^{-2}{B}^{-1/2}C{B}^{-1/2}\phi =W$$

(21a)

where the elements of matrix $C$ and vector $W$ are given by

$$C_{ij}=T_{x,ij}+\left({\displaystyle \frac{3}{8}}{\left({\displaystyle \frac{j^{\prime }\left(x_i\right)}{j\left(x_i\right)}}\right)}^2-{\displaystyle \frac{j^{″}\left(x_i\right)}{4j\left(x_i\right)}}+r_{\max }^2{j}^2\left(x_i\right)\left(U_{eff}\left(x_i\right)-E\right)\right){\delta }_{ij}$$

(21b)

(here, E is considered known) and

$$W_i=j^{1/2}\left(x_i\right)w\left(x_i\right).$$

(21c)

Furthermore, when solving DL equations, we have c = 1. For SE, this constant is determined by the wavefunction normalization condition

$${\displaystyle \underset{0}{\overset{r_{\max }}{\int }}{P}^{2}dr=1} \to c^2{r}_{\max }{\displaystyle \underset{0}{\overset{1}{\int }}{j}^2{y}^{2}dx}=c^2{r}_{\max }{\displaystyle \underset{0}{\overset{1}{\int }}{\varphi }^{2}dx}=1.$$

(22)

Finally, radial matrix elements involving a given function f(r) are evaluated as

$${\displaystyle \underset{0}{\overset{r_{\max }}{\int }}{P}_{k}f\left(r\right)P_{l}dr=}{c}_k{c}_l{r}_{\max }{\displaystyle \underset{0}{\overset{1}{\int }}{\varphi }_{k}f\left(x\right){\varphi }_{l}dx\approx }{c}_k{c}_l{r}_{\max }\mathrm{\Delta }x{\displaystyle \sum _{i=1}^N{\varphi }_k\left(x_i\right)f\left(x_i\right){\varphi }_l\left(x_i\right)}.$$

(23)

that is, they reduce to a simple summation over the N grid points.

3. Fisher Information Measure

The coordinate Fisher information measure is defined as [25,26,27,30]

$$I_r=4{\displaystyle \int \rho \left(r\right){\left|\overrightarrow{\nabla }\ln \rho \left(r\right)\right|}^{2}dV}=4{\displaystyle \int {\left|\overrightarrow{\nabla }\mathrm{\Psi }\left(r\right)\right|}^{2}dV}$$

(24)

where ρ(r) = |Ψ(r)|² stands for the spatial density distribution. Accordingly, the momentum space Fisher information is defined as,

$$I_p=4{\displaystyle \int \gamma \left(p\right){\left|\overrightarrow{\nabla }\ln \gamma \left(p\right)\right|}^2{d}^{3}p}$$

(25)

and γ is the corresponding momentum density distribution. Most of the present work deals with free Hydrogen atoms. Then, for a single particle system and a central potential, Equations (24) and (25) reduce to [40],

$$I_r=4\langle p^2\rangle -2\left(2\ell +1\right)\left|m\right|\langle r^{-2}\rangle$$

(26)

$$I_p=4\langle r^2\rangle -2\left(2\ell +1\right)\left|m\right|\langle p^{-2}\rangle$$

(27)

with m, the magnetic quantum number. For a unified treatment of all ℓ ≥ 0 levels, m is hereafter set equal to zero. Furthermore, <p²> is evaluated via the expression,

$${\langle p^2\rangle }_{n\ell }=2\left[E_{n\ell }-{\langle U\rangle }_{n\ell }\right].$$

(28)

Due to the absence of any confining potential in Equation (1), it is assumed that Equations (26)–(28) probe unconfined atoms. Of course, the imposition of a maximum radius (r_max) in the numerical solution of SE or DL is effectively, in its own right, equivalent to the system’s confinement within a box of this length. This has been a known issue for a long time in the community, and when the unconfined system is of interest, the choice of a sufficiently large value of r_max (to be practically unconfined) was always a matter of great concern and experimentation [1,2,21]. Computational limitations, on the other hand, dictate that r_max should not be unnecessarily large either. This is a central point of the present study.

Let us now consider the presence of a static electric field, where Equations (26) and (27) cannot be used because spherical symmetry no longer holds. In this case, our analysis shall be restricted to hydrogen’s 1s state (more generally, to non-degenerate s-states) and the coordinate Fisher information measure. We start from the system’s wavefunction, which, up to the first order correction, is written as,

$$\mathrm{\Psi }\left(r\right)=r^{-1}\left\{P_0^{\left(0\right)}{{\rm Y}}_0^0+F\cdot P_1^{\left(1\right)}{{\rm Y}}_1^0\right\}$$

(29)

where ${{\rm Y}}_{\ell }^m$ denote spherical harmonics, P_ℓ are radial wavefunctions, and F is the static electric field strength. For convenience, we set F = 1 throughout this work. The p-wave radial wavefunction $P_1^{\left(1\right)}$, obeys the following DL equation [21],

$$\left[-{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^2}{d{r}^2}}+U\left(r\right)+{\displaystyle \frac{1}{r^2}}-E^{\left(0\right)}\right]P_1^{\left(1\right)}=-3^{-1/2}r{P}_0^{\left(0\right)}$$

(30)

with E⁽⁰⁾, the unperturbed zeroth-order energy of the considered s-state. The static dipole polarizability a_d of that state is then related to its second-order energy correction [21],

$$a_d=-2E^{\left(2\right)}=-2\cdot 3^{-1/2}{\displaystyle \int P_1^{\left(1\right)}r{P}_0^{\left(0\right)}dr}.$$

(31)

Inserting Equation (29) into the definition (24) and taking into account (i) the DL Equation (30), (ii) the field-free SE obeyed by $P_0^{\left(0\right)}$ and (iii) the normalization of the latter, the space coordinate Fisher information measure is finally written as,

$$I_r=I_r^{\left(0\right)}+8\left\{{\displaystyle \int {\left(P_1^{\left(1\right)}\right)}^2\left[E^{\left(0\right)}-U\left(r\right)\right]dr}+{\displaystyle \frac{a_d}{2}}\right\}$$

(32)

where $I_r^{\left(0\right)}$ corresponds to the field-free s-state given by Equation (26). This is an interesting result, implying that the coordinate Fisher information measure and the static dipole polarizability are closely related.

In the following, our purpose is to examine the behavior of Equations (26) and (27) for static-field-free hydrogen atoms and Equation (32) under the static-field presence as a function of a number of computational parameters.

4. Results and Discussion

4.1. Schrödinger Equation: Fisher Information and Level Energy Systematics on r_max and N

We begin our analysis with a typical convergence procedure where energy eigenvalues are computed as a function of r_max, each time for fixed N. We concentrate on the more demanding ℓ = 0 states of the hydrogen atom and employ the transformation of Equation (11) in conjunction with the version of ${\widehat{T}}_x$ given by Equation (6). The exact hydrogenic bound state (negative) energy eigenvalues are given by E_n = −1/(2n²) a.u., where n is the principal quantum number. Over the total number of N eigenvalues and eigenvectors delivered by the solution of the eigenvalue equation, Equation (20d), the number of bound levels depends primarily on r_max and, to a lesser extent, on N (other eigenvalues correspond to continuum hydrogenic states and particle-in-a-box ones). Therefore, at present, we restrict ourselves to a state-by-state study. Figure 1 shows the typical evolution of the 1s ground state energy E_1s, coordinate, I_r, and momentum space, I_p, Fisher information measure as a function of r_max for the fixed value of N = 100. The exact values for these quantities are $E_{1s}^{\mathrm{exact}}$ = −0.5 a.u., $I_r^{\mathrm{exact}}$ = 4 a.u. and $I_p^{\mathrm{exact}}$ = 12 a.u. Note first that E_1s is higher than the exact value for small r_max (Figure 1a). As the latter increases, E_1s decreases continuously; it passes through the exact value and continues to decrease at a smaller rate. Consequently, the optimum r_max value for a given N cannot be determined variationally through the minimization (extremization in general) of the level energy itself. Variational minimization of the energy is evidenced if FGM is employed without any mapping (or, more precisely, with the linear mapping r(x) = r_maxx, 0 ≤ x ≤ 1). However, as expected from the preceding analysis, the emerging energy eigenvalues are quite inaccurate, and it is not worth discussing this case any further. Returning to the quadratic mapping of Equation (11), we observe that, contrary to the behavior of the energy as a function of r_max, there are optimal maximum radii, $r_{\max }^{\mathrm{opt}}$, for both I_r and I_p. The former quantity is minimized, its minimum value approaching its exact one (Figure 1b), while the latter is maximized, its maximum value approaching again its exact one (Figure 1c). The two extrema occur at slightly different $r_{\max }^{\mathrm{opt}}$ locations, and the minimum of I_r corresponds to a slightly more accurate energy eigenvalue. This observation persists for all ℓ ≥ 0 and n ≥ 1 levels. Although this point requires further investigation, a qualitative explanation can be based on the coordinate space probed by I_r and I_p, respectively. As is evident from Equation (28), I_r depends on <U> ∝ <r⁻¹>; that is, it is sensitive to the small-r quality of the radial wave function. On the other hand, I_p is proportional to <r²>, i.e., it is sensitive to the wavefunction’s large-r quality. The latter is obviously much more severely degraded for finite and small values of r_max.

By repeating the same type of calculation of Figure 1 for different numbers of grid points, we now associate in Figure 2 the computed Fisher information measure with the error in energy. It may be noted that the values of E_1s corresponding to the extrema of I_r and I_p (say, $E_{1\mathrm{s}}^{\mathrm{opt}}$) approach more and more $E_{1\mathrm{s}}^{\mathrm{exact}}$ as N increases. More importantly, the higher precision achieved is faithfully probed by these extrema; that is, smaller errors are associated with lower minima of I_r and higher maxima of I_p. Moreover, the inset of Figure 2 shows that the obtained $r_{\max }^{\mathrm{opt}}$ values increase (sub-linearly) with N and, therefore, optimum precision is obtained for specific ($r_{\max }^{\mathrm{opt}}$,N) pairs. In fact, it turns out that for a given precision on I_r, I_p and E_nℓ obtained for a chosen ($r_{\max }^{\mathrm{opt}}$,N) pair can be recovered by another ($r_{\max }^{\prime }$,N′) pair, such that N′ > N and $r_{\max }^{\mathrm{opt}}=r_{\max }^{\prime }{x}_N^2$, where x_N = (N − 1)/(N′ − 1). Furthermore, plots (a) and (b) of Figure 2 suggest error will be eliminated in the N→∞ limit and that the level’s optimal energy and maximum radius (the latter expected to be “infinite” or at least very large) may be found by an appropriate extrapolation to this limit (or the 1/N→0 one) [22,50]. Additionally, the fact that, for finite N, neither I_r nor I_p can probe a maximum radius, completely eliminating the absolute energy error, is to be attributed to the loss of mapped-FGM’s variational character. Nevertheless, the observed behavior of the Fisher information measure for the nonlinear mapping of Equation (11) seems to somehow restore this character, as long as the extremized quantity is not the level energy but, instead, the Fisher information measure itself.

Calculations show that when r_max < $r_{\max }^{\mathrm{opt}}$, and for either I_r or I_p, the computed energy eigenvalue sometimes fails to converge to the exact values in the N→∞ limit. More generally, as Figure 1 demonstrates, the relative error is higher when r_max < $r_{\max }^{\mathrm{opt}}$ with respect to the r_max > $r_{\max }^{\mathrm{opt}}$ case. Some of the causes of such a behavior may be elucidated through an inspection of Figure 3a, where E_1s is computed as a function of N for a fixed maximum radius r_max = 8 a.u. The latter value is sufficiently small for the hydrogen atom in this particular state to be considered loosely confined. Indeed, in this case, E_1s passes through the free atom’s $E_{1\mathrm{s}}^{\mathrm{exact}}$ value for N ≈ 150 and subsequently converges towards the exact energy eigenvalue of the confined atom [51] as N increases further. We may additionally note in Figure 3b,c that I_r and I_p exhibit extrema lying in the neighborhood of the free atom’s $E_{1\mathrm{s}}^{\mathrm{exact}}$. The conclusions that can be drawn from these observations are: (i) mapped-FGM always delivers the exact energy eigenvalues in the N→∞ limit and (ii) the extrema of I_r and I_p, as computed by Equations (26) and (27) and in conjunction with Equation (28), probe the energy eigenvalues of the unconfined system. These conclusions, after appropriate adaptations, might be useful in calculations of photoionization dynamics performed through the solution of the time-dependent Schrödinger equation (TDSE), where the spatial extension of a wave-packet varies with time and needs to be modified accordingly [52].

Let us now comment on the collective computation of several energy eigenvalues. We first note that the preceding state-by-state analysis verifies the expected increase in the corresponding $r_{\max }^{\mathrm{opt}}\left(n\right)$ values with increasing n. Moreover, as mentioned above, the number of bound levels emerging from such a calculation increases with r_max. Hence, one needs sufficiently large r_max values for computing many energy levels with acceptable precision, but, as Figure 1 suggests, these large r_max values will generally lead to a loss of accuracy for the ground and lower-lying levels. This behavior has long ago been predicted and observed [17,47]. On the other hand, conventional wisdom implies that, in any diagonalization procedure, the highest energy levels will exhibit the largest divergences from their exact values [41]. In fact, a cut-off (non-integer) effective principal quantum number v_c may be estimated based on the arguments presented earlier, i.e., by setting,

$$E_{\max }=-{\displaystyle \frac{1}{r_{\max }}}=-\frac{1}{2v_c^2} \to v_c={[r_{\max }/2]}^{1/2}.$$

(33)

The relative error of several members of the ℓ = 0 series is plotted in Figure 4a as a function of $n$, for fixed r_max = 136 a.u. and N = 100. Evidently, all our expectations are verified since the relative error monotonically decreases within the interval 1 ≤ n ≤ 6 and then gradually grows for n ≥ 7. The cut-off effective principle quantum number v_c ≈ 8.25 is also marked in Figure 4a and reasonably predicts the point of unacceptable loss of accuracy (i.e., n = 8, 9 and beyond, where the relative error deteriorates dramatically). It just misses the somewhat growing relative error of the n = 7 level.

Figure 4b shows the coordinate Fisher information difference $I_r-I_r^{\mathrm{exact}}$ corresponding to each $n$ = 1–8 level as a function of r_max. For the r_max` range shown, which includes the selected value of 136 a.u. for the computation of Figure 4a, minima exhibit the curves with n ≤ 6. In fact, this selected value corresponds to the optimum maximum radius for n = 6 and N = 100, where the $I_r-I_r^{\mathrm{exact}}$ curve has its minimum. Thus, by inspection of $I_r-I_r^{\mathrm{exact}}$ (or simply I_r(r_max)), one could predict the growing relative errors even for n = 7, and this criterion appears to perform somewhat better than that of Equation (33). An analogous behavior is observed for the maxima of I_p. As for the ℓ > 0 case, the same trend holds, although it appears that v_c does not depend solely on r_max (as Equation (33) implies), but it shows a small dependence on ℓ as well.

4.2. Dalgarno Lewis Equation: Fisher Information and Dipole Polarizability Systematics on r_max and N

The difference between the computed and the exact static dipole polarizability a_d of the ground state of hydrogen, as obtained by the mapped-FGM solution of Equation (30) and the integral of Equation (31), is plotted in Figure 5 as a function of r_max and for two values of N. Again, the mapping functions of Equation (11) and Equation (6) for ${\widehat{T}}_x$ are employed. The results show a behavior that is identical to what has emerged from the analysis of energy levels. Namely, a_d(r_max) does not show any extrema but rather simply different slopes before and after the exact value crossing. The coordinate Fisher information measure I_r, on the other hand, does show clear minima near the exact value.

As the number of grid points increases from N = 100 to N = 200, the minimum is deeper, and I_r converges towards the exact value (much like a_d). For N = 200 and at the minimum of I_r (i.e., for r_max = $r_{\max }^{\mathrm{opt}}$), the difference from the exact values for both I_r and a_d are of the order of 10⁻¹² a.u. This is comparable to the result given in [21], with the same mapping but for N = 500. The difference lies in the employed value of r_max that was selected via different arguments and not the extremization of the Fisher information measure. In fact, for N = 500 and a different radius with respect to the one chosen in [21] produces even more accurate results, the difference from the exact values being of the order of machine precision. Due precisely to this high accuracy obtained for both I_r and a_d, for a given specific s-state and sufficiently high N, the minimum of I_r becomes so shallow that it is occasionally difficult to locate, and then one simply searches for its stabilization within a pre-specified tolerance.

4.3. Schrödinger Equation: Fisher Information Systematics on Transformations h(x) = x^k, k ≥ 2

The preceding results demonstrated the ability of the Fisher information measure to point towards optimal ($r_{\max }^{\mathrm{opt}}$,N) pairs for the accurate calculation of energy levels and dipole polarizabilities. All calculations were performed with a single mapping function, h(x), as given by Equation (11). This mapping is clearly superior to the linear one, h(x) = x, as predicted by both the Fisher information measure and phase-space arguments [14] that were presented earlier. As a next step, let us now explore the possibility of the Fisher information measure to probe the relative computational efficiency offered by various mapping functions, possibly depending on one or several parameters. While both I_r and I_p may probe the optimum parameter values, our earlier results show that I_r performs slightly better. Therefore, for economy reasons, we focus our attention solely on I_r in this section. Furthermore, for properly confronting the Fisher information measure against phase-space notions, we choose to discuss the mapping function [12,17],

$$h\left(x\right)=x^k, k\ge 2,$$

(34)

for which, as it can be easily proved and the inset of Figure 6 shows, the ratio Λ ≡ A/A_box = π[2 r_max]^1/2/(2p_x,max) for the hydrogen atom continuously decreases as the integer k increases. Therefore, according to phase-space arguments, the efficiency of the mapping of Equation (34) is expected to degrade with increasing k. Earlier calculations, however, show the opposite trend, justifying the term “educated guess” that has been used in the description of this mapping [12,17].

For the 2s excited hydrogenic state, Figure 6 shows the coordinate Fisher information measure difference $\mathrm{\Delta }{I}_r\equiv I_r-I_r^{\mathrm{exact}}$ as a function of the energy difference $E_{2\mathrm{s}}^{\mathrm{exact}}-E_{2\mathrm{s}}$ for k = 2, 3, 4, 6 and 7. As in Figure 2, these differences are obtained by varying r_max. In fact, the same information can be obtained if ΔI_r is computed after replacing $I_r^{\mathrm{exact}}$ by I_r(r_max→∞) (in practice a sufficiently large radius r_max) and traced as a function of E_2s. For the shown k values, ΔI_r exhibits ever deeper minima with increasing k, thus reflecting the more accurate evaluation of E_2s. The important point made here is that I_r provides direct evidence of the performance enhancement afforded by the increase in k. More importantly, it allows for the demonstration that I_r is indeed able to probe this higher performance via its minimization, in stark contrast to parameter Λ. Similarly, I_p is also able to probe the superiority of k > 2 cases via its maximization. For s-states and all tested integer k values within the 2 ≤ k ≤ 7 interval, the only exception to the above behavior concerns the k = 5 case, where the role of I_r (and I_p) and E as a function of r_max is interchanged (i.e., E exhibits a minimum, but I_r does not). As mentioned earlier, this also holds for the linear k = 1 mapping. The reason for the discrepancy is not clear yet. It is, nevertheless, certain that the extremization of I_r and I_p is ensured for all even k values. Furthermore, even for k = 5, computational precision follows the trend and lies between k = 4 and k = 6. Some experimentation shows that for sufficiently large values of r_max, the absolute error |ΔX| = |X − X^exact| (X = E, I_r, I_p) is proportional to N^−k. Last but not least, as either N or k increases, the extrema of I_r and I_p are quite difficult to locate (with $I_r-I_r^{\mathrm{exact}}$ calculations reaching machine precision limit). In Figure 6, the formation of clear I_r minima is assisted by the very low (for an excited state) number of grid points N = 64. The above observations are particularly true for ℓ > 0 states, even for k = 2 and moderate N. Then, as mentioned earlier, one may search for the stabilization of I_r and I_p as a function of r_max within a pre-specified tolerance. On the other hand, this stabilization appears to be quite beneficial; since then, contrary to the situation depicted in Figure 4, the precision achieved for a large set of energy eigenvalues would be practically uniform for all levels until the energy cut-off, even for a necessarily large value of r_max.

Other mapping functions were also tested, with emphasis on those that assume the functional form x^k as x→0, while they are converted to x^k′, k′ < k for x > x_o, where k, k′ and the location x_o are treated as parameters to be optimized for given N. As an example among the tested mappings, we mention the function h(x) = x^k−1∙tan⁻¹(x/x_o)/tan⁻¹(1/x_o) (where k′ = k − 1). As it turned out and faithfully probed by the extrema exhibited by I_r and I_p in all cases, optimum results are achieved for x_o >> 1, where the mapping function tends towards the single k value mapping of Equation (34). As a last remark, a little experimentation shows that the mapping proposed by [14] and also adopted by [15], and which, in our case, writes h(x) = (x − atan⁻¹[bx])/(1 − atan⁻¹[b]), corresponds to Equation (34) with k = 3 when b = 1/a. With other combinations of parameters a and b, the mapping turns out to be either linear or difficult to use.

5. Conclusions and Outlook

We have presented a detailed investigation of the connections between the optimal mapped-FGM numerical solutions of hydrogenic Schrödinger and Dalgarno–Lewis equations and the extremization of coordinate and the momentum Fisher information measure. It has been demonstrated that, for a fixed number of grid points, the coordinate, I_r (momentum, I_p), Fisher information measure exhibits a minimum (maximum) as a function of the radial extension r_max of the mapped-FGM numerical solution. I_r probes the small-r quality of the radial wavefunction while I_p, the large-r one, while these two complementary quantities have been found to approach their exact values when higher precision is attested on the computed energy levels (or static dipole polarizabilities). The small investment of computational resources in this extremization procedure is amply “rewarded”, as Figure 1 proves, even in the simplest of problems, the hydrogen atom. Matrix diagonalization techniques for energy eigenvalues can benefit from the extra information and detect optimum parameters with more confidence. Additionally, relative error calculations in both Fisher measures succeeded in classifying mapped-FGM schemes (coordinate transformations) according to their efficiency and accuracy with a fixed number of grid points. The above arguments hold even when occasionally the extrema are hard to locate, as the Fisher information measure curves become practically flat at large r_max. Typical cases include sufficiently high numbers of grid points and ℓ > 0 states. Then, optimization can be simply and reliably based on the stabilization of the Fisher information measure, while this behavior can also be beneficial when many computed energy eigenvalues and, consequently, large radial extensions are required.

The framework laid out in the present work is, in principle, also applicable to Alkali metals or any other atoms with a single valence electron when excitation of only this valence electron is of concern. Then, the many-body problem is reduced to a single-active-electron one by describing the ionic core system by a local radial model potential. The valence electron is subjected to this potential, which exhibits a Coulomb tale at large distances, while it is modified at short distances to simulate the core’s inner electrons. Promising preliminary results in analytic, parametric, generally ℓ-dependent, radial model potentials (such as the one provided in [53] for the Li atom), as well as the Hulthén potential (for ℓ = 0 [54]), suggest that the present methodology (validated here in the hydrogen atom) can, indeed, be successfully applied to these more complicated atomic systems, whose systematic investigation will be the subject of future work.

When free (unconfined) atoms are of interest, the maximum radius of the mapped-FGM calculation is a rather restrictive parameter since there is always the possibility of affecting numerical results through confinement. The behavior of both Fisher information measures provides a very effective way to circumvent the problem. Yet, it would be quite beneficial for the evolution of FGM to follow earlier suggestions where the semi-infinite radial interval is mapped into a finite one [23,24]. Then, parametrization as a function of the radial extension would be eliminated, and the Fisher information measure could be used to optimize the values of a number of other parameters included in the mapping function. Going one step beyond, one could set as the “ultimate goal” for a future extension of the present work the development of a more general extremization procedure based on the extreme physical information principle (see, for example, [32,55,56,57]) in which the potential energy curve and any related information would be the only required input and a potential-specific mapping function would be the output. We are currently working toward this goal.