How to numerically solve a 1D time-independent Schrödinger equation for two interacting particles

Question

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?

Answer 1

As @yarchik pointed out, this is strictly speaking not a two-electron problem, but a problem of "distinguishable particles," since the antisymmetry postulate (and spin) is ignored. Moreover, the 1D Coulomb potential in the initial equations can't be just $1/x$ -- it has to have an absolute sign because otherwise it would be attractive on one side and repulsive on the other. In the actual code, we're dealing with a different potential which has the correct symmetries, and also doesn't have the pathological divergence at the origin.

However, you have to keep in mind that both particle coordinates can be positive or negative, and therefore the boundaries of the NDEigensystem calculation cannot be chosen at $x_1=x_2=0$.

Instead, you have to choose the boundaries symmetrically. Of course the goal should be to choose the boundary distance d large enough so that the Dirichlet conditions on that boundary don't affect the results. This makes the calculation quite time-consuming. With the settings (d=10) of the first attempt in the question, this gives the following results:

{egv1,egs1}=With[{d=10,n=3},NDEigensystem[{1.5 f[x1,x2] + V2[x1,x2] 
  f[x1,x2] - 1/2 Laplacian[f[x1,x2],{x1,x2}],
   DirichletCondition[f[x1,x2]==0,True]},f,{x1,-d,d},{x2,-d,d},
    n,Method->{"SpatialDiscretization"->{"FiniteElement",
     {"MeshOptions"->{"MaxCellMeasure"->0.001}}},"Eigensystem"->{"Arnoldi"}}
]];  

With the choice d=20, the result for the lowest-energy state shouldn't look very different because the images above indicate that the solution has indeed decayed almost to zero, so that it won't feel the Dirichlet condition very much. However, the computation ran too long for me to wait for the result. For the excited states, there seems to still be a significant effect of the boundary, so it may be necessary to increase the simulation domain.