Solution of the 3D logarithmic Schrödinger equation with a central potential

Some form of the time-independent logarithmic Schrödinger equation (log SE) arises in almost every branch of physics. Nevertheless, little progress has been made in obtaining analytical or numerical solutions due to the nonlinearity of the logarithmic term in the Hamiltonian. Even for a central potential, the Hamiltonian does not commute with L 2 or L z ; the Hamiltonian is invariant under the parity operation only if the wave function is an eigenstate of the parity operator. We show that the solutions with well-defined parity can be expressed as a linear combination of eigenstates of L 2 and L z , where the parity restrictions on l and m determine the nodal structure of the wave function. The dominant contribution in the sum is designated as l ˜ and m ˜ ; these serve as approximate quantum numbers. Using an iterative finite element approach, we also carry out fully converged numerical calculations in 1D, 2D and 3D for the special case of a Coulomb potential. The nodal structure of the wave functions and the expectation values L 2 and L z 2 are consistent with the analytical predictions. Values for the energy and expectations values are tabulated for the low-lying states. The methods developed for this problem are applicable across many areas of physics.


Introduction
The time-independent logarithmic Schrödinger equation (log SE) [1][2][3] is a non-linear modification of the standard Schrödinger equation, where y S r ; ( )  is the wave function and > S 0 is a constant parameter. An excellent review of the origins and applications of the logarithmic wave equation is given in [4]. The log term is associated with entropy and ensures a factorization of wave functions for composed systems [1]. It is distinct from the potential and cannot do dynamical work on the system [5]. (Other forms of the non-linear SE also include a source term that depends on the probability density of the wave function; see e.g., [6].) The log SE equation is notoriously difficult to solve. Standard variational methods fail because the Hamiltonian is non-Hermitian. The logarithmic term cannot be treated with time-independent perturbation theory because it dramatically changes the structure of the ordinary Schrödinger wave function.
The log SE appears in many fields of physics, including quantum optics [7], nuclear physics [8], transport and diffusion phenomena [9], dispersion studies [10], stochastic quantum mechanics [11], Bose-Einstein condensation [5,12], quantum gravity [12][13][14] and superfluids [15,16]. For example, it was recently shown that the logarithmic term y S r ln ; | ( )|  of equation (1) dominates the standard Gross-Pitaevskii terms for Helium-4 at low temperatures in certain regimes [16]. A superfluid vacuum theory using a log SE appears in one formulation for the Higgs potential [17,18]. The log SE even has geophysical applications in magma transport [19]. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
We recently considered the special case of the log SE with a Coulomb potential and arbitrary S [20]. We obtained an exact analytical solution for the ground state. We developed an iterative finite element approach to solve the 1D log SE for states with spherical symmetry (where = l 0 and = m 0 are good quantum numbers). Herein, we extend our work in two directions. In section 2, we analyze the symmetries and nodal structures of the bound states for an arbitrary central potential, In section 3 we discuss the iterative finite element method and the additional challenges of extending it to higher dimension. In section 4, we solve the log SE for the Coulomb potential in 1D, 2D, and 3D. Concluding remarks are given in section 5.

Analytic properties of the solutions to the log SE for a central potential
The logarithmic Schrodinger equation can be written in spherical coordinates q r, and f, We consider the case where V r ( ) is an attractive potential. (The wave function is dependent on the parameter S, but we suppress the index hereafter.) Atomic units are used throughout. Because of the presence of the logarithmic term, the Hamiltonian no longer commutes with L 2 or L .
z The log SE for the central potential is not separable in any coordinate system. The wave functions are not orthogonal but obey where n and n′ represent the collective quantum numbers of two different physical states. (See appendix A of [20] for a derivation.) The state y n is orthogonal (in the normal sense) to y k which are solutions of The states y k , which satisfy y y d á ñ = n k nk , | , do not exist physically and they are not solutions of equation (2) unless = k n. However, they do exist mathematically and will appear in any iterative solution of the log SE. (This is discussed further in the next section.) In the absence of the logarithmic term, the Hamiltonian is invariant under the parity operation r r,   and the wave functions are eigenstates of the parity operator, With the addition of the logarithmic term, the Hamiltonian remains invariant under the parity operation if and only if the wave functions satisfy equation (5). For the states with well-defined parity, the eigenfunctions of L 2 and L z still represent an appropriate basis set to expand the wave function. That is, This suggests that the (real) bound state wave function can be expressed as a linear combination of products of associated Legendre polynomials q P cos n the summations, l is even(odd) if l˜is even(odd); m is even(odd) if m is even(odd). The approximate quantum numbers l˜and m refer to the dominant contribution to the sum. The third quantum number n is associated with radial excitations. Thus, the traditional designation ¼ s s p p p 1 , 2 , 2 , 2 , 2 x y z ( ) serves as an approximate classification scheme. The contributing states for   l 0 3 (and the predicted nodal surfaces) are summarized in table 1. Note that nd z 2 is unique in that the state = l 0 contributes to the sum; this is the only state where < l l . This has important consequences on the value of á ñ L 2 and the structure of the wave function; there are no nodal surfaces corresponding to a fixed value of θ and the wave function does not vanish at the origin. Similar examples appear for  l 3. Although equation (7) is useful for understanding the nodal structure of the wave functions, it is impractical from a computational point of view. One cannot solve a set of coupled 1D equations for the radial functions R r nl m lm ( )˜because the coupling term between different l m , ( ) angular momentum states involves the logarithm of the absolute value of the wave function, which is itself a sum of angular momentum states weighted by the radial function. In other words, one must solve equation [2] directly. We will confirm by rigorous numerical calculations that all of the numerical solutions obey the parity condition of equation (5), and the predictions for the nodal structures given in table 1 are accurate.

Iterative finite element method for non-linear Hamiltonians
In order to solve the log SE for a central potential, we use the iterative finite element approach that was introduced in [20] for the 1D log SE. There is a vast literature on the finite element method (FEM) [21], so we only review the important features of the iterative scheme. In all of the results reported here, we use a uniform discretization of the grid; this allows for a systematic convergence study. For purposes of this discussion, the single number n is used to label the energy and wave function associated with a particular level. The log SE (for a given value of S) must be solved for each value of n separately, because y n appears in the Hamiltonian itself. In the first iteration, the asymptotic solution is used in the logarithmic term; the unknown wave function y q f r, , where δ is the tolerance on the iterative process. In most cases, this simple approach begins to diverge after a few iterations.
There are two sources of instability. First, a slight change in the location of a nodal structure has a dramatic effect on the contribution of the logarithmic term. The matrix element y y y á ñ + +  If the location of a nodal structure changes too dramatically in a single iteration, the energy may begin to oscillate about the exact value with increasing amplitude. In order to enhance stability, the wave function (after a few initial iterations) is adjusted slowly by using a weighted average of the last two iterations. In the ]/ in the log term. The linear combination is renormalized, although it is already very close to unity. By increasing N, the convergence rate is slowed but the   m l m 2 , even; f m cos( ) iterative process is stabilized. This process does not affect or restrict the symmetry of the final wave function. In the 2D and 3D calculations, it was often necessary to carry out hundreds of iterations.
The second source of instability is due to the presence of the unphysical states, which are solutions to equation (4). The desired state y n can go through one or more avoided crossings with adjacent unphysical states. The steps described above can help prevent a pre-mature mixing of the levels, but there is no way to eliminate the effect of these states completely. When the states begin to mix, the iterative process will begin to diverge. To monitor the early onset of this mixing, we evaluate the overlap integral y y á ñ

|˜( )
If the overlap integral begins to decrease, the iterative process must be stopped immediately. The accuracy of the energy is limited by this point of 'closest approach' to the avoided crossing.

Results for the Coulomb potential
For purposes of illustrating the method, we consider the special case of the attractive Coulomb potential This causes a compression of the wave function in the radial direction, even for the highly excited states. Excitations in the radial direction will raise the energy significantly.
In order to proceed in a systematic way and to check the accuracy in higher dimensions, we solve equation (8) in 1D (spherical symmetry, = L 0), 2D (azimuthal symmetry, = m 0), and 3D. In all cases, we do not impose a priori any additional parity or symmetry requirement on the wave functions. We only impose the asymptotic boundary condition, The value of r max is systematically increased until the energy is stable to the desired number of digits. In the 3D calculations, we also require that the wave function (and its derivatives) be singled valued at f f p = = 0 and 2 .

1D solution of the log SE for spherical symmetry
In this special case, = l 0 and = m 0 are good quantum numbers. The ground state (which we label s 1 ) is exactly solvable [20], and the solution is The lowest three sstates are obtained numerically with the iterative FEM. The number of 1D elements is n r and the size of each element is r n .

2D solution of the log SE for azimuthal symmetry
For the special case where the states are independent of the azimuthal angle, = m 0 is still a good quantum number. We solve the Coulomb log SE on a 2D grid in r and q. The number of elements is´q n n r and the area of each element is ṕ q r n n . / Except for s 1 and p 2 z states, the accuracy of the final answer (for a given grid) is limited by the onset of an avoided crossing with an unphysical state. It is important to test the convergence by varying q n even for the s states as this impacts the energy of the unphysical states involved in the avoided crossings. The energy of the 2s state lies above that of the p p d 2 , 3 and 3 z z z 2 states. This is due to the energy cost of a radial excitation. In 2D, the s 2 state has an avoided crossing with the unphysical state below it, which has a p 2 state like structure. Thus, even using the same radial grid as the 1D calculation, the accuracy of the s 2 energy is less because the iterative process is stopped when the overlap integral begins to decrease. The d 3 z 2 and p 3 z states are particularly problematic, as an avoided crossing occurs after the first iteration with the level below it; the solution oscillates between states of two symmetries and the overlap y y á ñ  Table 3.   increases with increasing n. This accounts for the dramatic increase in á ñ L 2 from p 2 z to p 3 . z For the 3d z 2 state, the non-zero value for Ã n2 0˜h as two important consequences. First, á ñ L 2 for d 3 z 2 is lower than the value for the p 3 z state even though the = l 2 contribution is dominant. In addition, the wave function is non-zero at the origin; the nodal curve is a function of both r and q, as predicted in section 2.

3D solution for the log SE with a Coulomb potential
In the 3D calculations, the number of elements is´q   previous calculations, the s 1 state is free of stability problems; there are no nodal surfaces or avoided crossings with unphysical states. The logarithmic term in the Hamiltonian removes the degeneracy between the s 2 and p 2 z state. However, the degeneracy between the p p 2 , 2 x y and p 2 z states remains. For a fixed grid, the FEM energy for p 2 z is slightly lower than p 2 x and p 2 y because the wave function is essentially exact in the azimuthal coordinate; any numerical error due to the approximation of the wave function raises the energy. The p 2 z state converges with no problems. The next two levels, p 2 x and p 2 y , are also numerically degenerate. A linear combination of these states is not a solution of the log SE because the Hamiltonian is non-linear. The iterative process immediately begins to oscillate between the two states; that is, if y n i corresponds to p 2 y , then y + n i 1 corresponds p 2 .
x Once again, it is necessary to jump across the avoided crossing early in the iterative process.
The results for the 3D iterative FEM calculations of the log SE are given in table 4. For the three p 2 states, the expectation values á ñ á ñ x y , 2 2 and á ñ z 2 are interchanged exactly as expected. The p 2 energy levels are converged to five digits; we record three digits for the expectation values for the purpose of comparison between the three states. These results confirm that the spatial symmetry of the hydrogen p 2 states is not destroyed by the logarithmic term; the approximate quantum numbers are an appropriate scheme for labeling the energy levels and the corresponding wave functions. All of the 3D solutions are eigenfunctions of the parity operator; i.e., they obey equation (5).

Conclusion
For a central potential, we have shown that the wave functions for the logarithmic Schrodinger equation can be expanded as a sum over angular momentum states; the dominant contribution is designated l  and m; these play the role of approximate quantum numbers. Requiring that the wave functions are eigenstates of the parity operator restricts the allowed states in the expansion; this in turn determines the nodal structures in the wave  functions. The analysis for a central potential will be useful for studying Bose-Einstein condensates in spherically symmetric atom traps [22][23][24][25] as well as spherically symmetric nuclear models [8].
We solved the log SE for a Coulomb potential using the iterative FEM. In spite of the logarithmic term, the bound state spectrum is very much like the atomic hydrogen solutions. This is a completely unexpected result. The exact quantum numbers for the hydrogen atom assume the role of approximate quantum numbers, as predicted by our analysis for any central potential. For the case of zero angular momentum, it was previously found that linear combinations of soliton solutions known as Gaussons [20] could accurately model the radial functions of excited states for the Coulomb potential. Herein, we have shown that for non-zero angular momentum, the wave functions can be expressed as a linear combination of spherical harmonics and radial functions, based on parity arguments. If these radial functions can also be modelled with Gaussons, this suggests that any resolution of the atomic system could be feasibly extended to the molecular problem using the tools of computational quantum chemistry, e.g., the LCAO (Linear Combination of Atomic Orbitals) approach. The importance of a nonlinear superposition property for the log SE is already under investigation [26]. Solving a logarithmic version of the Schrödinger-Poisson equation could even be contemplated to further cosmological investigations into dark matter [27,28].
In conclusion, we anticipate that the iterative FEM developed for the solution of the log SE will have widespread applications in many fields of physics. The approach is easily adapted to different coordinate systems and there is no restriction on the form of the physical potential.