NDSolve aborted

I'm trying to simulate the formation of vortices in a BEC condensate with Mathematica. This problem was solved a long time ago (see this paper). After making several experiments and changes in the equations and boundary conditions I believe that I need a larger spatial resolution to "see" the vortices. Using a low resolution only gives some resemblance of vortices. However, when I try to increase the resolution either Mathematica aborts the calculation or the memory grows without control. I'd tried the method described in this question without success.

The full code is:

clearStatus[] := showStatus[""];
clearStatus[];
a = 20;
xl = yl = -a;
xr = yr = a;
tmax = 1.;
α = 0.01;
ϵ = 0.08;
ClearAll[vortex];

vortex[x_, y_] := (1/Sqrt[6 π ] ) Exp[-(x^2 + y^2)/(6)];

eqn = I (1 + α I)*
Derivative[0, 0, 1][ψ][x,y,t] == -0.5 Laplacian[ψ[x, y, t], {x, y}] +
I 0.8 (x D[ψ[x, y, t], {y,1}] - y D[ψ[x, y, t], {x,1}])
+ (600 Abs[ψ[x, y, t]]^2 +
(x^2 + (1+ϵ) y^2)/2)*ψ[x, y, t];

bcs = {ψ[xl, y, t] == vortex[xl, y], ψ[xr, y, t] == vortex[xr, y],
ψ[x, yl, t] == vortex[x, yl], ψ[x, yr, t] == vortex[x, yr]};

(* bcs={ψ[xl,y,t] == ψ[xr,y,t],ψ[x,yl,t] == ψ[x,yr,t]};*)

ics = ψ[x, y, 0] == vortex[x, y];

nxy = 200;
sol = First[NDSolve[{eqn, ics, bcs}, ψ, {x, xl, xr}, {y, yl, yr},
{t,0,tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> nxy, "MaxPoints" -> 20 nxy,
"DifferenceOrder" -> "Pseudospectral"}},
EvaluationMonitor :> showStatus["t = " <> ToString[CForm[t]]]]];


I want to solve exactly this equation, with possible different or zero parameters α and ϵ (dissipation and anisotropy). The result must contain vortices, which must appear as holes in "the pie" shown below. Any help will be appreciated. Thanks.

Observation: Please note that the answer given below is not the physical solution to my problem, so I emphasize to preserve the differential equation of this post.

Why do not you use equation (5.11) from the article cited? In this case, vorticity appears quickly enough and exists long enough. See the code and examples

    a = 4;
xl = yl = -a;
xr = yr = a;
tmax = 2;
\[CapitalOmega] = 1.7; \[Beta] = 100; \[Gamma] = -2;
vortex[x_, y_] := ((x + I*y)/Sqrt[\[Pi]]) Exp[-(x^2 + y^2)/2];
eqn = Derivative[0, 0, 1][\[Psi]][x, y, t] ==
0.5 Laplacian[\[Psi][x, y, t], {x, y}] -
I \[CapitalOmega] (x*D[\[Psi][x, y, t], {y, 1}] -
y*D[\[Psi][x, y, t], {x,
1}]) - (\[Beta] Abs[\[Psi][x, y, t]]^2 + \[Gamma]*(
x^2 + y^2))*\[Psi][x, y, t];
bcs = {\[Psi][xl, y, t] == vortex[xl, y], \[Psi][xr, y, t] ==
vortex[xr, y], \[Psi][x, yl, t] ==
vortex[x, yl], \[Psi][x, yr, t] == vortex[x, yr]};
ics = \[Psi][x, y, 0] == vortex[x, y];
sol = NDSolveValue[{eqn, ics, bcs}, \[Psi][x, y, t], {x, xl, xr}, {y,
yl, yr}, {t, 0, tmax}]; // Quiet
{Grid[{{"\[CapitalOmega]=", \[CapitalOmega]}, {"\[Beta]=", \[Beta]}, \
{"\[Gamma]=", \[Gamma]}}],
Table[ContourPlot[Abs[sol], {x, xl, xr}, {y, yl, yr},
PlotLabel -> Row[{"t=", t}], PlotRange -> All, Contours -> 20,
PlotPoints -> 50], {t, 0, tmax, .2*tmax}]}


I tested the method for NDSolve[], which allows investigating the solution of this problem in a wide range of parameters. Here is an example of the formation of quantum vorticity in an electromagnet field:

Lx = 4;
Ly = 4;

tmax = 2;
\[CapitalOmega] = 1.3; \[Beta] = 100; \[Gamma][t_] :=
Sin[2*Pi*t]/Sqrt[2];
vortex[x_, y_] := Exp[-(x^2 + y^2)/2];
eqn = Derivative[0, 0, 1][\[Psi]][x, y, t] ==
0.5 Laplacian[\[Psi][x, y, t], {x, y}] -
I \[CapitalOmega] (x*D[\[Psi][x, y, t], {y, 1}] -
y*D[\[Psi][x, y, t], {x,
1}]) - (\[Beta] Abs[\[Psi][x, y, t]]^2 - \[Gamma][
t - x - y]^2)*\[Psi][x, y, t];
bcs = {\[Psi][-Lx, y, t] == vortex[-Lx, y], \[Psi][Lx, y, t] ==
vortex[Lx, y], \[Psi][x, -Ly, t] ==
vortex[x, -Ly], \[Psi][x, Ly, t] == vortex[x, Ly]};
ics = \[Psi][x, y, 0] == vortex[x, y];
sol = NDSolveValue[{eqn, ics, bcs}, \[Psi][x, y, t], {x, -Lx,
Lx}, {y, -Ly, Ly}, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 40, "MaxPoints" -> 100,
"DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6];
Table[ContourPlot[Abs[sol], {x, -Lx, Lx}, {y, -Ly, Ly}, Mesh -> None,
Contours -> 20, PlotLegends -> Automatic,
PlotLabel -> Row[{"t=", t}], ColorFunction -> Hue,
AspectRatio -> Automatic, PlotRange -> All,
FrameLabel -> {"x", "y"}], {t, 0, tmax, .2*tmax}]


• Thank you for your solution, the sign in the potential do not agree with my code, I'll check it. It seems in the paper the potential as the same sign as the non-linear term as well. However, your solution looks correct. – Gluoncito Jun 28 '18 at 12:38
• The authors consider the harmonic attracting potential, whereas I checked the harmonic repulsive potential due to the rotation of the liquid. I could not reproduce their examples. However, my examples are similar to them. – Alex Trounev Jun 28 '18 at 14:14
• In my first comment, I didn't realize that the eq 5.11 of the cited paper is not the OP equation. It is the equation for the "normalized gradient flow" method. I'm glad that the vortices shows up, but I'm frustrated for not being able to obtain them from a direct simulation (from eq. 2.5). I suppose Mathematica must be able to perform the integration and that I'm missing something else. – Gluoncito Jun 29 '18 at 0:11
• The rotation is expressed by the angular momentum operator, the harmonic potential only serves to confine the condensate. In order to perform the simulatation with the equation you used it is necesary to renormalize the wave function in each timestep (and change the sign of the potential). This is why we don't get the same solutions. – Gluoncito Jun 29 '18 at 0:19
• If you want to compare your calculations with the results in Fig. 6.1 you should use (5.11). See 6.1. Initial data for computing ground state. … Here we present a 2D example to evolve the CNGF (5.11) with its BEFD discretization for four diﬀerent initial data. – Alex Trounev Jun 29 '18 at 3:02