The Schrödinger–Newton Equation
Introduction
The Schrödinger equation predicts that the wave packet for any object will invariably spread out over time, leading to a delocalized center of mass. This dispersion is in conflict with the everyday experience of well-localized classical objects, whose center of mass always seems to have a well-defined position. The Schrödinger–Newton (SN) equation offers a possible solution to this predicament by including the effects of gravitational self-interaction. As proposed by Lajos Diósi [1], this approach introduces a self-generated potential that counteracts the natural tendency of the wave packet to spread. For objects with sufficient mass, this self-gravity would be enough to prevent dispersion, allowing stationary and localized solutions [2]. Furthermore, the SN equation has gained significant interest in recent years, as it represents the consistent non-relativistic limit of semi-classical gravity, directly connecting to the foundational question of whether gravity must be quantized [3]. Investigating these stationary solutions is therefore crucial for exploring gravity's potential role in localization and for analyzing the behavior of this fundamental semi-classical framework.
This application example will be based on the work by Bernstein et al. [4], who build upon the work by Jones [5]. In the latter, the self-generated gravitational potential was proposed as a universal "localizing" term, and the SN equation is derived in the context of a many-body problem. For simplicity, the focus will be on a cold Bose–Einstein condensate (BEC) composed of 
identical bosons. This physical system allows for a crucial simplification known as the Hartree ansatz, where the total N-particle wavefunction is approximated as a product of identical single-particle wavefunctions, as shown in Equation (1).
Furthermore, the analysis is restricted to computing S-wave eigenstates, which are spherically symmetric. This symmetry is mathematically advantageous, as it reduces the three-dimensional problem to a one-dimensional one, thereby making the numerical solution considerably more manageable.
With this ansatz, it is possible to show that the SN system is the one shown in Equations (2) and (3), for the eigenfunction 
and potential 
, where 
[4].
Furthermore, by choosing appropriate units and scaling factors, it is possible to rewrite the SN system as in Equations (4) and (5).
Then, by making the standard transformation 
, the equations to solve in spherical coordinates become:
Note that Equation (6) is a Poisson-type equation. However, it will be referred to as the Newton equation here, since it represents Poisson’s equation for Newtonian gravity and helps distinguish it from the Schrödinger–Poisson equation that makes use of a Poisson-type equation used in electrostatics. In contrast, the Schrödinger–Poisson arises when the goal is to find a wavefunction consistent with the electrostatic potential energy generated by itself.
Going Step by Step
The strategy used to solve the problem is iterative. While it is possible to combine Equations (7) and (8) into a single integro-differential equation, this results in a nonlinear equation that is nontrivial to solve. Dividing the problem into two equations allows one to solve only one equation at a time and then use the result as input for the other. Repeating this process until the eigenvalue 
converges leads to the correct eigenfunction 
and potential 
.
The geometry under consideration is very simple: 
ranges from 
to an external radius 
. The external radius must be large enough for the boundary conditions 
and 
to be physically meaningful—that is, for the functions 
and 
to have effectively decayed to zero at that distance.
For these reasons, the steps to be taken are the following:
1. Set the outer boundary radius 
.
2. Supply an initial guess for 
.
3. Solve the Newton equation (9) for 
.
4. Solve the Schrödinger equation (10) for the 
th eigenvalue and eigenfunction.
5. Iterate the previous two steps until the eigenvalue converges sufficiently.
6. Iterate the previous five steps, increasing 
until 
stops changing appreciably.
To illustrate this process, steps 1 through 4 will be shown one by one first, and later, helper functions will be given to approach steps 5 and 6.
Step 1: Set the Outer Radius 
An initial value for 
is guessed, and its value will be revised later. Also, a mesh for the region is defined in accordance.
Needs["NDSolve`FEM`"]
ℛ = 100;ToGradedMesh is used to generate the mesh because it allows a higher element density in the region where the wavefunction—and likewise the gravitational potential—attain their largest values, rather than in the region where both exhibit a decaying behavior.
mesh = ToGradedMesh[{Line[{{0}, {ℛ}}], <|"Alignment" -> "Left", "ElementCount" -> 100|>}]Set Up the Operators
Then the Schrödinger operator that specifies the eigenvalue problem and the left-hand side of the Newton equation are defined.
ℋ = SchrodingerPDEComponent[{u[r], {r}}, <|"ReducedPlanckConstant" -> 1, "SchrodingerPotential" -> ϕ[r]|>]nOp = Laplacian[ϕ[r], {r, θ, ϕ}, "Spherical"]Step 2: Supply an Initial Guess for 
u0[r_, a_] := a r Exp[-(r^2/10)]
uini = Part[u0[r, a] /. NSolve[Integrate[ u0[r, a]^2, {r, 0, ℛ}] == 1, {a}], 2]Is important to remember the substitution 
was made, so 
.
Plot[{(uini/Sqrt[4π]r), uini}, {r, 0, 30}, PlotRange -> All, PlotLegends -> {"f", "u"}]NIntegrate[uini^2, {r, 0, ℛ}]Step 3: Solve the Newton's Equation for 
Now Equation (11) is solved, where the initial wavefunction guess is specified on the right-hand side of the equation.
Φ = NDSolveValue[{nOp == (uini^2/r^2), DirichletCondition[ϕ[r] == 0, r == ℛ]}, ϕ, {r}∈mesh]Step 4: Solve the Schrödinger Equation for the 
th Eigenstate
Using the Newtonian potential obtained in the previous step, the Schrödinger eigensystem is then solved. The "Shift" option for the "Arnoldi" method is used, and a reordering of the eigenstates is required to ensure that the lowest eigenvalue appears first.
sortStates[s_] := SortBy[s, Re]{en, funs} = sortStates@NDEigensystem[{ℋ /. ϕ -> Φ, u[0] == 0, u[ℛ] == 0}, u[r], {r}∈mesh, 5, Method -> {"Eigensystem" -> {"Arnoldi", "MaxIterations" -> 3000, "Shift" -> -10, "BasisSize" -> 50}}]For illustrative purposes, the focus here is on the ground state, but any eigenvalue could be analyzed.
Plot[(funs[[1]]/Sqrt[4 π]r), {r, 0, 30}, PlotRange -> All]Note that the overall sign of the wavefunction is physically irrelevant, since only its modulus squared carries physical meaning.
Step 5: Iterate the Previous Two Steps until the Eigenvalue Converges
For this next step, it is necessary to solve the same equations over and over. For that reason, helper functions are defined.
First, a helper function is defined to solve the Schrödinger equation. It takes the gravitational potential as an input, along with the total number of eigenvalues requested and the mesh to be used. Having the mesh as an input for the helper function will be useful for step 6 of the process, in which the external radius 
is gradually increased, requiring the mesh to be updated accordingly.
solveEigenvalue[Φ_, nT_, mesh_] := sortStates@NDEigensystem[{ℋ /. ϕ -> Φ, DirichletCondition[u[r] == 0, True]}, u[r], {r}∈mesh, nT, Method -> {"Eigensystem" -> {"Arnoldi", "MaxIterations" -> 3000, "BasisSize" -> 80, "Shift" -> -10}, "Interpolation" -> {"ExtrapolationHandler" -> {Automatic, "WarningMessage" -> False}}}
]Wolfram Language functions like NDEigensystem typically select the optimal numerical method automatically. However, for certain problems, you may need to adjust method options to improve the solution's quality or, in this case, the eigenvalue's convergence.
NDEigensystem operates by first discretizing the PDE and then using the Eigensystem function to solve the result. The options for Eigensystem that can be used within NDEigensystem are detailed in this link to the documentation. For the problem at hand, increasing the "BasisSize" is specially helpful for step 6, in which the external radius 
is progressively increased.
The "Interpolation" option is specified to allow extrapolation. This enables the function to evaluate the solution at points 
, which is required for step 6, as will be clear later.
In a similar way, a helper function to solve Newton's equation is defined.
solveNewton[u_, mesh_] := NDSolveValue[{nOp == (u^2/r^2), DirichletCondition[ ϕ[r] == 0, r == mesh["Bounds"][[1, 2]]]}, ϕ, {r}∈mesh]Then a pair of helper functions are defined to perform steps 3 and 4 until the change in the energy eigenvalues is less than one part in 
. They take as an input the external radius 
, the eigenstate of interest n, the total number of eigenvalues to calculate nt> n, and the initial guess for the wavefunction 
.
schrodingerNewtonStep[previousStep_List, mesh_, n_, nT_] := Module[{ϵ, u, ϕ, ϵOld, uOld, ϕOld},
{ϵOld, uOld, ϕOld} = previousStep;
ϕ = solveNewton[uOld, mesh];
{ϵ, u} = solveEigenvalue[ϕ, nT, mesh][[{1, 2}, n]];
{ϵ, u, ϕ}
]schrodingerNewtonSolver[ℛ_, n_, nT_, uini_] := Module[{energies, u, ϕ, mesh, steps, tol = 10^-9},
mesh = ToGradedMesh[{Line[{{0}, {ℛ}}], <|"Alignment" -> "Left", "ElementCount" -> 8000|>}];
steps = NestWhileList[
schrodingerNewtonStep[#, mesh, n, nT]&, {Infinity, uini, None}, Abs[#1[[1]] - #2[[1]]] > tol&, 2];
{energies, u, ϕ} = Transpose[steps];
{Rest[energies], Last[u], Last[ϕ]}
]Note that the element count for the mesh has been increased. This is because, as will be clear later, when the radius 
is increased, more elements are needed.
AbsoluteTiming[{enList, 𝓊, Φ} = schrodingerNewtonSolver[ℛ, 1, 5, uini]]GraphicsRow[ListPlot[#]& /@ {enList, Differences@enList}]
Plot[(𝓊/Sqrt[4 π]r), {r, 0, 30}, Rule[...]]Plot[Φ[r], {r, 0, 30}, ...]Step 6: Iterate the Previous Five Steps, Increasing 
until 
Stops Changing
The equation is then solved for a larger radius, 
, using as the initial guess the solution obtained for the smaller 
. This is the reason the extrapolation in the solution 
was allowed earlier.
schrodingerNewtonSolver[ℛ + 1, 1, 5, 𝓊]A slightly higher eigenvalue is obtained by increasing the radius from 
to 
.
enList//LastNow find the outer radius 
for which the eigenvalue 
is stable.
ℛ0 = ℛA new helper function is defined to compute the energy successively, as in the previous step, but for progressively larger outer radii 
. The process continues until the difference between consecutive eigenvalues is smaller than one part in 
.
nextRadiusStep[previousStep_, n_, nT_] := Module[{nextℛ, ϵ, u, ϕ, ϵOld, uOld, ϕOld, ℛOld},
{ϵOld, uOld, ϕOld, ℛOld} = previousStep;
{ϵ, u, ϕ} = schrodingerNewtonSolver[ℛOld, n, nT, uOld];
nextℛ = ℛOld + 10 Log[1. + 10^8Min[1, Abs[Last[ϵ] - ϵOld]]];
{Last[ϵ], u, ϕ, nextℛ}
]
findConvergedBoundary[n_, nT_] := Module[{tol = 10^-6, ℛ = ℛ0, ℛOut, steps, energies, u, ϕ},
{energies, u, ϕ, ℛOut} = NestWhile[ nextRadiusStep[#, n, nT]&, {Infinity, uini, None, ℛ}, Abs[#1[[1]] - #2[[1]]] > tol&, 2]
]
{timeTaken, {ϵn, 𝓊, Φ, ℛOut}} = AbsoluteTiming[ findConvergedBoundary[1, 5]]timeTaken / 60.The ground state eigenvalue agrees with the values reported by Bernstein et al. [4] and Harrison et al. [6], namely 
. Moreover, as illustrated below, the plots of the ground state eigenfunction and the corresponding gravitational potential closely match those presented by Bernstein et al. [4].
Plot[(𝓊/Sqrt[4 π]r), {r, 0, 30}, ...]Plot[Φ[r], {r, 0, 30}, ...]Conclusions
This analysis has successfully computed stationary solutions to the SN equation, and the methodology shown here can be used to analyze higher states instead of only the ground state. The types of eigenstates found here can be used as initial conditions for studies of the time evolution and stability of the time-dependent version of the SN system. One example of what is possible is the work by Harrison et al. They mathematically tested the stability of these states by calculating how they respond to tiny disturbances, finding that only the ground state is stable. Subsequently, they used these same stationary solutions as the starting point in a full numerical simulation to observe how the unstable higher states decay and settle into the stable ground state over time [6].
Additionally, it has been shown here how using high-level functions like NDSolveValue and NDEigensystem makes it straightforward to build the complex, iterative procedures required for this kind of analysis.
References
[1]: Diósi, L. (1984). "Gravitation and Quantum-Mechanical Localization of Macro-Objects." Physics Letters A, 105(4–5), 199–202. https://doi.org/10.1016/0375-9601(84)90397-9.
[2]: Salzman P J (2005) "Investigation of the Time Dependent Schrödinger-Newton Equation." PhD Dissertation, the University of California at Davis. http://www.dirac.org/physics/dissertation.
[3]: Bahrami, M., Großardt, A., Donadi, S., & Bassi, A. (2014). "The Schrödinger–Newton Equation and Its Foundations. New Journal of Physics, 16(11), 115007. https://doi.org/10.1088/1367-2630/16/11/115007.
[4]: Bernstein, D. H., Giladi, E., & Jones, K. R. W. (1998). "Eigenstates of the Gravitational Schrödinger Equation." Modern Physics Letters A, 13(29), 2327–2336. https://doi.org/10.1142/S0217732398002473.
[5]: Jones, K. R. W. (1995). "Gravitational Self-Energy as the Litmus of Reality. Modern Physics Letters A, 10(08), 657–667. https://doi.org/10.1142/S0217732395000703.
[6]: Harrison, R., Moroz, I., & Tod, K. P. (2003). "A Numerical Study of the Schrödinger Newton Equations." Nonlinearity, 16(1), 101–122. https://doi.org/10.1088/0951-7715/16/1/307.