Solving the 1D single-electron time-independent Schrödinger equation has been demonstrated using NDEigensystem
here. There, the single-electron Schrödinger equation is
$$ \left(-\frac{1}{2}\frac{d^2}{dx^2}-\frac{1}{|x|}\right)\psi (x) =E\psi(x) $$
So is it possible to use NDEigensystem
to find a few of the lowest eigenstates of a two-particles Schrödinger equation
$$
\left(-\frac{1}{2}\frac{d^2}{dx_1^2}-\frac{1}{2}\frac{d^2}{dx_2^2}-\frac{1}{|x_1|}-\frac{1}{|x_2|}+\frac{1}{\left|x_1-x_2\right|}\right)\psi (x_1,x_2)=E \psi(x_1,x_2)
$$
Since we are ignoring the spin here, the two particles are distinguishable. To make things simpler, I replace the Coulomb potential with the soft-core one:
$$ \left(-\frac{1}{2}\frac{d^2}{dx_1^2}-\frac{1}{2}\frac{d^2}{dx_2^2}-\frac{1}{\sqrt{1+x_1^2}}-\frac{1}{\sqrt{1+x_2^2}}+\frac{1}{\sqrt{1+(x_1-x_2)^2}}\right)\psi (x_1,x_2)=E \psi(x_1,x_2) $$
V[x1_, x2_] := -(1/Sqrt[1 + x1^2]) - 1/Sqrt[1 + x2^2] + 1/Sqrt[1 + (x1 - x2)^2]
With[{d = 10},Plot3D[V[x1, x2], {x1, -d, d}, {x2, -d, d}]]
Then solve the eigenstates the same way as that in the one-electron problem:
{egv1, egs1} =
With[{d = 10,n=3},
NDEigensystem[{1.5 f[x1, x2] + V[x1, x2] f[x1, x2] -
1/2 Laplacian[f[x1, x2], {x1, x2}],
DirichletCondition[f[x1, x2] == 0, True]},
f, {x1, 0, d}, {x2, 0, d}, n,
Method -> {"SpatialDiscretization" -> {"FiniteElement",
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}},
"Eigensystem" -> {"Arnoldi"}}]]; // AbsoluteTiming
(* {41.5109, Null} *)
With[{d = 10},
Plot3D[#[x, y], {x, 0, d}, {y, 0, d}, ImageSize -> Medium,
PlotRange -> All] & /@ egs1]
However, when I change the solving range from d = 10
to d = 20
for a larger range, the solutions change
{egv2, egs2} =
With[{d = 20,n=3},
NDEigensystem[{1.5 f[x1, x2] + V[x1, x2] f[x1, x2] -
1/2 Laplacian[f[x1, x2], {x1, x2}],
DirichletCondition[f[x1, x2] == 0, True]},
f, {x1, 0, d}, {x2, 0, d}, n,
Method -> {"SpatialDiscretization" -> {"FiniteElement",
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}},
"Eigensystem" -> {"Arnoldi"}}]]; // AbsoluteTiming
(* {410.335, Null} *)
Plot3D[#[x, y], {x, 0, d}, {y, 0, d}, ImageSize -> Medium,
PlotRange -> All] & /@ egs2
Compare the first solution for these two cases we have
Show[{Plot3D[egs1[[1]][x, y], {x, 0, 10}, {y, 0, 10},
ImageSize -> Medium, PlotRange -> All, PlotStyle -> Red],
Plot3D[egs2[[1]][x, y], {x, 0, 20}, {y, 0, 20}, ImageSize -> Medium,
PlotRange -> All, PlotStyle -> Green]}]
We can see that the first eigenstates in these two cases peak at different positions in space. Their eigenvalues are also different:
egv1
egv2
(* {1.4056, 1.42235, 1.60257} *)
(* {1.26956, 1.26974, 1.33544} *)
This seems to mean the original solutions are not correct, since we are trying to find the bound states of the system which should not depend on the box size.
So how can I get the correct solutions?