Solution of the logarithmic Schrödinger equation with a Coulomb potential

The nonlinear logarithmic Schrödinger equation (log SE) appears in many branches of fundamental physics, ranging from macroscopic superfluids to quantum gravity. We consider here a model problem, in which the log SE includes an attractive Coulomb interaction. We derive an analytical solution for the ground state energy and wave function as a function of the strength of the logarithmic interaction. We develop an iterative finite element method to solve the Coulombic log SE for the spherically symmetric states. The ground state results agree with the exact solution to better than one part in 1010. The excited states (n>1) are converged to better than one part in 108. We also construct a remarkably simple variational wave function, consisting of a sum of Gaussons with n free parameters. One can obtain an approximation to the energy and wave function that is in good agreement with the finite element results. Although the Coulomb problem is interesting in its own right, the iterative finite element method and the variational Gausson basis approach can be applied to any central force Hamiltonian.

Bialynicki-Birula and co-workers investigated the log SE in the context of nonlinear wave mechanics [11][12][13]. Recently there has been more widespread interest in this highly nonlinear form of the Schrödinger equation. It has been hypothesized that superfluid vacuum theory (SVT) might be responsible for the mass mechanism (in contradistinction to the Higgs boson or perhaps working in tandem with it). Others have proposed connections between quantum gravity and Bose-Einstein condensates in the form of quantum liquids [14]. A version of SVT favors a logarithmic Schrödinger equation over a Gross-Pitaevski equation [9,15]. In particular, the logarithmic model can account for the upside-down Mexican hat shape of the Higgs potential and its solution is claimed to be even more stable and energetically favorable than the model with a quartic (Higgslike) potential [9,16,17]. Although the Higgs boson reported at 125 Gev has been confirmed, the Higgs potential has not been firmly established.
The reported resemblance between the gravitational dilaton and the Higgs boson [18] becomes all the more compelling because the log SE is a common feature. Thus, in the present work, we consider precise methods, both analytical and computational, for accurately solving the log SE. This not only serves as a sequel to our previous work in dilatonic quantum gravity [19] but provides the tools by which to make a number of alternate and future investigations with the log SE.
Herein, we consider a computational model consisting of a Coulomb potential with an additional logarithmic term for the following reasons: (i) The electromagnetic interaction appears in proposed models for the Higgs potential [17].
(ii) The Newtonian and post-Newtonian gravitational potentials in the general relativistic model have the form 1/r n [20,21]; n=1 corresponds to the case considered here.
(iii) The two-body dilatonic quantum gravity model [22], consisting of only one spatial dimension, leads to Coulombic terms for the more general dilatonic problem of [19] in three spatial dimensions.
The outline of this work is as follows. First, we present the log SE governing our model and the unusual nonorthogonality condition that results from the logarithmic term. The log SE with a Coulomb interaction has a simple analytic solution only for the ground state. This exact solution provides us with a benchmark by which to test our numerical methods. Second, we present the computational details of the iterative FEM approach [23][24][25]. The finite element method has been successfully used to solve the Schrödinger equation for a variety of atomic and molecular systems, in 1D, 2D and even 3D. Based on a variational principle, it provides rigorous upper bounds to the energy. In particular, the FEM method succeeded in isolating the effect of a singularity structure at the Fermi level (10 −15 m range) within a continuous structure at the atomic level (10 −10 m range) [26]. The FEM solutions for the log SE with the Coulomb interaction are presented for the ground state and a number of excited states.
We also present variational results using linear combinations of (time-independent) Gaussons [27], that is, products of exponentials and Gaussians. (These should not be confused with the time-dependent Gaussons that are the free-particle solitonic solutions of the log SE.) Despite the simplicity of the basis, the energies and wave functions are remarkably close to the fully converged FEM solutions. Concluding remarks, commentary and discussion are presented at the end.

Theoretical model
The time-independent logarithmic Schrödinger equation with a Coulomb interaction can be expressed (in atomic units) as ; , , , ) is a constant parameter, and S e( ) is the total bound state energy. The system supports bound states only for S0. The log SE is not separable in the radial and θ coordinates. However, since e ln 0 im = f | | , m remains a good quantum. We consider here the spherically symmetric case, where the wave function is independent of the angular coordinates, ψ (S; r, θ , f )=R(S; r) Y 0 0 (θ , f where n is the principle quantum number and E S S ln 4 n n S 2 The Hamiltonian operator is nonlinear because it depends on the absolute value of the radial function, R S r ; n | ( )|. Note that the difference between E n and ε n is a direct consequence of the 3D nature of the problem. Experimentally, it is differences between energy states that are measured and for the same value of S, E n and ε n differ only by a constant (offset). This is true only for the special case when the wave function is independent of the angular coordinates. Moreover, if we make the rescaling R r E Sf r exp n = -( ) ( ) ( ), the energy E n (S) is eliminated [12, equation (18)] from (2)(although E n appears in the normalization of f (r)). Hereafter, we refer to E n as the energy, while the total bound state energy of the system is The (non-degenerate) radial eigenfunctions of equation (2) are not orthogonal [11]. As shown in appendix A, they obey: where the integration is over the radial coordinate. The logarithmic interaction is pathologically non-linear: , and singular at the nodes of the radial function.
Using the Maple computer algebra system [28], the leading asymptotic solution for equation (5) assuring square integrability for both the solution and its first derivative is given by: where C is an arbitrary real constant. Note that R S r R S r ; ; n n = | ( )| ( ) for r 0 < < ¥. Thus, asymptotically R(r) has the general form: where b may be a function of S. The functional form of equation (7) is that of a Gausson, i.e. the product of a Gaussian in r of the form r exp ). Note, however, that since the Gaussian dominates the asymptotic behavior even for small S, the coefficient b need not be positive. Equation (7) expresses the asymptotic form of all the bound eigensolutions.

Ground state
Equation (2) is analytically solvable only for the ground state, where where erfc is the complementary error function. Note that R 1 (S; r) of equation (9) is consistent with the asymptotic form of equation (7) for a single Gausson with b=1. In the limit S E N 0, 1 2, 2 S 1  =-= and R 1 =exp (−r) which is indeed the ground state of the hydrogen atom. Using equation (8), and setting 0 , we find through symbolic manipulation with Maple that the ground state energy E 1 reaches a maximum when S is governed by: . The integer k refers to the particular branch. In this case, it is not the main branch at k=0 but rather k=−1 for which W(−1, x) is real valued on the interval [−e −1 ,0]. As predicted earlier [19,22,31], the Lambert W function does make its appearance in this body of work. Solving equation (11) leads to: S E 0.640 250059733 0.302 028005290 12 The relevance of this maximum (which also occurs for excited states) is not yet clear, since the bound state energy ε n has no such feature. Nevertheless, it reflects a 'competition' between the Coulomb potential and the logarithmic term.

FEM computational aspects
In order to gain insight into the analytical structure of the excited states, we solved equation (2) iteratively using the finite element approach. Details of the finite element method are given in [23][24][25], so we provide only a brief summary here. In 1D finite element analysis, the continuum is truncated at r max , where the function and its derivative are required to vanish (we suppress the index n and the parameter S for clarity). The continuum 0<r<r max is discretized into nel small regions called elements. In each element iel, the unknown radial wave function R(r) is locally approximated by a sum of six 5th degree polynomials, The polynomials f j (x) have the property that the six unknown expansion coefficients c j iel are the value of R(r) and dR/dr at the two endpoints and the midpoint of the iel element. Using equation (13) in equation (2) and projecting onto the basis functions, one obtains a 6×6 local matrix equation for each of the nel elements, With the exception of the logarithmic term, all of the integrals are simple polynomials and can be evaluated exactly. We used a 16-point Gauss quadrature to evaluate the logarithmic term. (We also carried out calculations at higher quadrature, to verify that the numerical integration was converged to the desired accuracy; see appendix B.) Imposing continuity of the wave function and its derivative across element boundaries (i.e., requiring that c c ) results in a generalized eigenvalue problem, where the two matrices are of order 4×nel with a bandwidth of 11. The resultant generalized eigenvalue problem was solved for the n th eigenvalue only using a banded storage LAPACK routine DSBGVX. The solution yields the energy E n and the expansion coefficients c j iel from which one can reconstruct a piecewise continuous function R n (S; r). For a given n, we find only the n th eigenvalue and eigenfunction; the lower eigenvalues l=1, 2, K, n−1 are unphysical and correspond to the states which satisfy Note that R n appears in the log term, not R l . These unphysical states obey the normal orthogonality condition since the effective potential is the same. There are only two parameters in the finite element calculation; the value of r max and the number of elements. This enables one to study convergence in a very systematic way. The FEM calculation is converged if increasing r max (and keeping the element size the same) or increasing nel (and keeping r max the same) does not change the value of the energy. The results of a typical convergence study are given in appendix B.
For the ground state calculation, the logarithmic term is ignored in the zeroth order iteration; that is, the equation is solved for the pure hydrogenic case. In the i th iteration i>0, the solution R i 1 1 is used in the logarithmic term. The energy is considered converged at the i th iteration when E E 10 For the excited states, this simple iterative scheme is not sufficient due to the pathological nature of the logarithmic term at the nodes. A slight change in the position of the node moves the singularity. Although the convergence is initially monotonic, the energy begins to oscillate about the exact solution, and in some cases, the oscillations increase in amplitude. To correct this problem we use a weighted average of the last two iterations. For the zeroth iteration, the logarithmic term is ignored. For 0<i10, the solution R n i 1 is used in the logarithmic term (as described above for the ground state). For i>10, the weighted average (I−i) R i−2 /I+i R i−1 /I is used in the logarithmic term. The value of I is increased until the convergence is monotonic. In general, the maximum number of iterations needed for convergence is roughly half the value of I. We emphasize that the value of I does not impact the final value for the energy, only the rate of convergence. For n>1, the energy is considered converged at the i th iteration when E where N(S) is a normalization constant, b 1 (S) and b 2 (S) are free parameters, and d 1 (S) is uniquely determined by the nonorthogonality rule of equation (4). Although equation (16) is not the exact solution, it provides an excellent approximation to the FEM wave function over the full range r 0 < < ¥. We further note that an error in the wave function of order δ gives rise to an error in the energy of the same order δ. This is in contrast to the linear Hamiltonian case, where the error in the energy is δ 2 . We note that the proof for the linear operator relies on the orthogonality of the wave function, which does not hold in this case.
In calculating the expectation value of the Hamiltonian R R 2 2  á ñ | | , all of the integrals can be evaluated exactly with the exception of the logarithmic term, which leads to integrals of the form r e d e dr ln 1 , 17 where k is a positive integer and B is some combination of b 1 and b 2 . For the n=2 excited state, d 1 is negative and thus the integral in (17) is a Cauchy principal value. The nonorthogonality condition between the n=2 state and the ground state also leads to such integrals. The class of integrals of (17) have no known analytical solutions and the derivatives with respect to the free parameters are even more unwieldy.
As an alternative, we use the hybrid symbolic-methods [38][39][40], in particular the generation of optimized MATLAB code [41] to apply the variational principle. We used MATLAB functions to numerically solve all the integrals of the class in (17) but all other integrals are solved analytically using Maple. Since derivatives of the integrals in (17) with respect to the parameters b 2 and d 1 are problematic, to find the minimum energy, we compute the expectation values for the energy R R 2 2  á ñ | | over a fine grid of values b 1 and b 2 with d 1 being determined from the nonorthogonality condition of (4). We then select the minimum, subject to the constraint that N ensures that R 2 is normalized to unity.
The nested loops over the parameters b 1 and b 2 are computationally expensive. However, the double loop over grid values of b 1 and b 2 can be reduced to a single loop over b 1 , by using the value of the n=2 node obtained from the FEM results. Requiring  (4). This leaves only n free parameters, b i , i=1, 2, K , n; however, the nested loops over the parameters b i become computationally prohibitive for n>2. We therefore employ the method described above for n=2, and use the known value of the (n−1) nodes to collapse the problem to a single loop over b 1

Finite element results
Our goal is to obtain energy levels that are accurate to better than one part in 10 10 for n=1, and better than one part in 10 8 for n>1. All calculations were carried out in double precision. In table 1, we compare the FEM results with the exact energy E 1 (S) given in equation (8). FEM results for n=2, 3 and 4 are given in tables 2-4, respectively; we include both the energies E n (S) and the position of the nodes. There is a well-defined maximum in E n (S) , (but not in ε n (S)). The energies for n=1, 2, 3, and 4 are plotted as a function of S in figure 1. One important check on the accuracy of the FEM wave functions was the nonorthogonality condition. We stress that we did not enforce equation (4) a prior, but found each state in a separate calculation. In table 5, we compare the value of R r R r n m á ñ ( )| ( ) with equation (4). It is an important feature of this approach that the FEM wave functions automatically satisfy the nonorthogonality condition.
In figure 2, we show the radial function for n=4 for S=0, 0.5, 1, 2 and 5. The logarithmic term is responsible for a dramatic compression of the wave function towards the origin.

Variational principle
The results of the variational calculation, i.e. the parameters for the ansatz of (16), for the n=2 state are tabulated in table 6 as a function of S; the values are rounded to 6 digits. These were obtained according to the method in subsection 3.2 using the nodes from the FEM results as tabulated in table 2. Note that these results   were confirmed without using the nodes by a more laborious calculation over grids for both the b 1 and b 2 parameters essentially confirming the same outcome. As S gets larger in magnitude, absolute discrepancies between the FEM and variational results increase, although the relative error decreases: for S=5, Δ E=0.02 (0.6 %) and for S=0.05, Δ E=10 −4 (2 % ). As we can see from table 6, the b 1 coefficient becomes negative at around S≈0.05 and d 1 becomes smaller in magnitude as S approaches zero, thus ensuring the correct S 0  limit. Figure 3 shows how well the variational function fits the FEM wave function to within plotting accuracy for three different values of S, i.e. S=1, 2 and 5. In figure 4, we plot the parameters as function of S. The magnitude of b 1 and b 2 appear to grow almost logarithmically with S In figures 5 and 6, we compare the variational wavefunction obtained from the ansatz equation (18) with the FEM wave function for the states n=3 and n=4. We used the nodes given in tables 3 and 4 to reduce the number of free parameters to one. It is clear that the variational wave function is within plotting accuracy of the accurate FEM result. As expected, the variational energies lies slightly above the FEM energies.

Discussion
In summary, we have obtained an analytical solution to the ground state for the log SE with a Coulomb interaction. For the spherically symmetric case, the FEM method yields extremely accurate solutions for the lowlying bound states. A simple variational wave function consisting of a sum of n Gaussons is shown to provide an excellent approximation for the excited states. This is particularly important because the most developed industry for quantum chemistry calculations employs Gaussians and Slater type functions (see [38][39][40] and   references therein). This suggests that many of the computational tools developed in quantum chemistry may be relevant even in the presence of the highly non-linear logarithmic term. We have considered the spherically symmetric case with a Coulomb potential. From here, a number of investigations with our numerical methods can be contemplated. First, the central potential used in the FEM calculation can easily be modified; for example, one could consider potentials that arise in quantum gravity or the Higgs potential, as discussed in the Introduction. Since the asymptotic behavior is dominated by the logarithmic term, it is expected that the Gaussons will continue to serve as an excellent variational basis. An equally challenging problem is to consider the 1D continuum states of the log SE with a Coulomb interaction (or other central potential).
In order to go beyond the spherically symmetric case, 2D FEM can be used to solve the log SE in the variables r and θ for fixed m, which remains a good quantum number for a central potential. In the 2D case, the logarithmic term diverges along the nodal lines. The number of nodal lines will increase with increasing energy. In the limit S 0  , the nodal lines must approach the correct hydrogenic limit, and the two new quantum numbers must collapse into n and ℓ.
Finally, the iterative approach developed here, which ensures monotonic convergence, could potentially be applied to the Gross-Pitaevski equation or to other nonlinear differential eigenvalue equations in field theory.  - Taking the inner product of (A.1) with R m and (A.2) with R n , one obtains , and exploiting the Hermiticity of the kinetic and potential energy operators, only two terms remain: Since E m >E n , the nonorthogonality condition is given by =´=for S=1.
The only parameters in the FEM calculation are the value of the cutoff radius r max and the size of the elements. The convergence study for the ground state at S=1 is shown in table B1. We also verify the accuracy of the Gauss quadrature. With only four elements, the energy is accurate to one part in 10 6 . For the excited states, it is necessary to introduce a weighted average of the last two iterations in order to obtain monotonic convergence. The weighting parameter does not affect the final value of the energy, only the rate of convergence. Results for the state n=3 and S=1 are shown in table B2. The convergence with respect to the Gauss quadrature is not strictly monotonic because the accuracy is slightly dependent on location of the sampling points relative to the nodes.