# Solving Stochastic Gross-Pitaevskii equation

I am trying to solve the Stochastic Gross-Pitaevskii equation from this paper https://arxiv.org/pdf/1409.0146.pdf. But I have no idea how to solve adding a noise term. I like to see the wave function after adding a temperature similar to https://arxiv.org/pdf/2107.06680.pdf fig(1). Below is the problem statement from https://arxiv.org/pdf/1409.0146.pdf and the results I want to see.

ClearAll["Global*"]

xtab = Table[i/10, {i, 1, 300}];
ttab = Table[i/10, {i, 1, 301}];

Noisetab = Table[Random[Real, {-1/5, 1/5}], {301}];

\[Eta] = Interpolation[
Join[Table[{-xtab[[i]], ttab[[i]], Noisetab[[i]]}, {i, 300}],
Table[{0, ttab[[i]], Noisetab[[i]]}, {i, 300}],
Table[{xtab[[i]], ttab[[i]], Noisetab[[i]]}, {i, 300}]]];

hbar = 1; mass = 1;
\[Gamma] = 0.01; g = 0.1; temp = 13.89; \[Mu] = 22.41; \[Omega] = 0.5;
\[Epsilon] = 1 - I*\[Gamma]; kb = 1;

pse[x_, t_] := Exp[-0.2 x^2];
sol = NDSolveValue[{ -
D[u[x, t], t] == -0.5 \[Epsilon] D[u[x, t], {x, 2}] +
0.5 \[Epsilon] \[Omega]^2 x^2 u[x,
t] + \[Epsilon] g Abs[u[x, t]]^2 u[x, t] - \[Epsilon] \[Mu] u[
x, t] + Sqrt[2  \[Gamma] kb temp] \[Eta][x, t], u[x, 0] == 0,
u[30, t] == 0, u[-30, t] == 0},
u, {x, -30, 30}, {t, 0, 20},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 161, "MaxPoints" -> 201,
"DifferenceOrder" -> 4}}, MaxSteps -> 10^6];

DensityPlot[Abs[sol[x, t]]^2, {t, 0, 20}, {x, -30, 30},
AspectRatio -> 1/2, Frame -> True, FrameTicks -> Automatic,
PlotPoints -> 200, ImageSize -> 1000,
LabelStyle -> {24, Bold, Large, Black},
ColorFunction -> "BlueGreenYellow",
FrameLabel -> {{Style["x", FontFamily -> "Times New Roman",
FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30],
None}, {Style["t", FontFamily -> "Times New Roman",
FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30],
None}}]

I tried pseudospectral Method

ClearAll["Global*"]

xtab = Table[i/10, {i, 1, 800}];
ttab = Table[i/10, {i, 1, 801}];

Noisetab = Table[Random[Real, {-1/5, 1/5}], {801}];

\[Eta] = Interpolation[
Join[Table[{-xtab[[i]], ttab[[i]], Noisetab[[i]]}, {i, 800}],
Table[{0, ttab[[i]], Noisetab[[i]]}, {i, 800}],
Table[{xtab[[i]], ttab[[i]], Noisetab[[i]]}, {i, 800}]]];

hbar = 1; mass = 1;
\[Gamma] = 0.01; g = 0.1; temp = 13.89; \[Mu] = 22.41; \[Omega] = 0.5;
\[Epsilon] = 1 - I*\[Gamma]; kb = 1;

pse[x_, t_] := Exp[-0.2 x^2];
sol = NDSolveValue[{
I D[u[x, t], t] == -0.5 \[Epsilon] D[u[x, t], {x, 2}] +
0.5 \[Epsilon] \[Omega]^2 x^2 u[x,
t] + \[Epsilon] g Abs[u[x, t]]^2 u[x, t] - \[Epsilon] \[Mu] u[
x, t] + Sqrt[2  \[Gamma] kb temp] \[Eta][x, t], u[x, 0] == 0,
u[30, t] == 0, u[-30, t] == 0},
u, {x, -30, 30}, {t, 0, 80},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 100, "MaxPoints" -> 161,
"DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6];

DensityPlot[Abs[sol[x, t]]^2, {t, 0, 80}, {x, -30, 30},
AspectRatio -> 1/2, Frame -> True, FrameTicks -> Automatic,
PlotPoints -> 200, ImageSize -> 1000,
LabelStyle -> {24, Bold, Large, Black},
ColorFunction -> "BlueGreenYellow",
FrameLabel -> {{Style["x", FontFamily -> "Times New Roman",
FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30],
None}, {Style["t", FontFamily -> "Times New Roman",
FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30],
None}}]

The result looks like

Can anyone help me with phase coherence g(x_1,x_2) and its dynamics? Also how to get a ground-state solution like https://arxiv.org/pdf/2107.06680.pdf with the inclusion of temperature? I tried with 2D equation from Solving logarithmic 2D Nonlinear GPE Equation.

ClearAll["Global`*"]
\[Gamma] = 0.01; g = 0.1; temp = 13.89; \[Mu] = 22.41; \[Omega] = 0.02;
\[Epsilon] = 1 - I*\[Gamma]; kb = 1;

data = Flatten[
Table[{{x, y, t}, Random[Real, {-1/100, 1/100}]}, {x, -40,
40}, {y, -40, 40}, {t, 0, 10}], 2];
\[Eta] = Interpolation[data, InterpolationOrder -> 1];
boundary = 40; xl = yl = -boundary; xr = yr = boundary;
finalt = 10;
seedwave[x_, y_] :=
Exp[-0.05 (x^2 + y^2) +
I*svalue ArcTan[x, y]] (x^2 + y^2)^(svalue/2);
pnumber =
NIntegrate[
Exp[-0.05 (x^2 +
y^2)], {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \
\[Infinity]}];(*particle number check*)
svalue = 0; lvalue = 0;

sol = NDSolveValue[{-D[\[Psi][x, y, t],
t] == -0.5 \[Epsilon] Laplacian[\[Psi][x, y, t], {x,
y}] + \[Epsilon] Abs[\[Psi][x, y, t]]^2 \[Psi][x, y, t] Log[
Abs[\[Psi][x, y, t]]^2] - \[Epsilon] \[Mu] \[Psi][x, y, t] +
Sqrt[2  \[Gamma] kb temp] \[Eta][x, y, t], \[Psi][xl, y, t] ==
0, \[Psi][xr, y, t] == 0, \[Psi][x, yl, t] ==
0, \[Psi][x, yr, t] == 0, \[Psi][x, y, 0] == 0.1}, \[Psi][x, y,
t], {x, xl, xr}, {y, yl, yr}, {t, 0, finalt},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 40, "MaxPoints" -> 81,
"DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6]
Plot3D[Evaluate[Abs[sol]^2 /. t -> finalt], {x, xl, xr}, {y, yl, yr},
PlotRange -> All, MeshStyle -> White, ColorFunction -> "Rainbow",
PlotTheme -> "Marketing", PlotPoints -> 100]

And the phase

Plot3D[Evaluate[Arg[sol] /. t -> finalt], {x, xl, xr}, {y, yl, yr},
PlotRange -> All, MeshStyle -> White, ColorFunction -> "Rainbow",
PlotTheme -> "Marketing", PlotPoints -> 100]

• What's wrong with the error message you get? That needs to be fixed first. Nov 3, 2022 at 10:28
• Due to the stochastic term, you might need to use ItoProcess instead of NDSolve, which would entail manually discretizing your equation in space. Here's an example applied to a stochastic reaction-diffusion equation. Nov 3, 2022 at 13:28
• Thanks for the reference @Chris K. Nov 3, 2022 at 14:05
• @ArghaDebnath I see you've edited your question to use ItoProcess now. Did it work or are you still having problems? Nov 4, 2022 at 13:02
• No, @Chris K it is running for a long time with no answer. Nov 4, 2022 at 16:08

We can compute correlation function with OrnsteinUhlenbeckProcess[] using mean in time as follows

tmax = 60; xmax = 60; nmax = 240;
pR = RandomFunction[
OrnsteinUhlenbeckProcess[0, 1., xmax/nmax/Sqrt[2]], {0, xmax,
xmax/nmax}]; pI =
RandomFunction[
OrnsteinUhlenbeckProcess[0, 1., xmax/nmax/Sqrt[2]], {0, xmax,
xmax/nmax}];
xR = Interpolation[pR, InterpolationOrder -> 1]; xI =
Interpolation[pI, InterpolationOrder -> 1];
eta[x_] = xR[x] + I xI[x];
tR = RandomFunction[
OrnsteinUhlenbeckProcess[0, 1., tmax/nmax/Sqrt[2]], {0, tmax,
tmax/nmax}]; tI =
RandomFunction[
OrnsteinUhlenbeckProcess[0, 1., tmax/nmax/Sqrt[2]], {0, tmax,
tmax/nmax}];
etR = Interpolation[tR, InterpolationOrder -> 1]; etI =
Interpolation[tI, InterpolationOrder -> 1];
eta1[t_] = etR[t] + I etI[t];

Solution

hbar = 1; mass = 1;
\[Gamma] = 0.01; g = 0.1; temp = 13.89; \[Mu] = 22.41; \[Omega] = \
0.5;
\[Epsilon] = 1 - I*\[Gamma]; kb = 1;
sol = NDSolveValue[{I D[u[x, t],
t] == -0.5 \[Epsilon] D[u[x, t], {x, 2}] +
0.5 \[Epsilon] \[Omega]^2 x^2 u[x,
t] + \[Epsilon] g Abs[u[x, t]]^2 u[x,
t] - \[Epsilon] \[Mu] u[x, t] +
Sqrt[2 \[Gamma] kb temp] eta[x + xmax/2] eta1[t]/4,
u[x, 0] == 0, u[xmax/2, t] == 0, u[-xmax/2, t] == 0},
u, {x, -xmax/2, xmax/2}, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 141, "MaxPoints" -> nmax,
"DifferenceOrder" -> 4}}, MaxSteps -> 10^6] // Quiet;

Visualization

DensityPlot[Abs[sol[x, t]]^2, {t, 0, tmax}, {x, -20, 20},
AspectRatio -> Automatic, PlotPoints -> 100,
ColorFunction -> "LightTemperatureMap", FrameLabel -> {"t", "x"},
PlotRange -> All, PlotLegends -> Automatic]

Correlation function

phi[x_] = Sum[sol[x, t], {t, tmax - 20, tmax, 1/5}]/100;
int[x_, x1_] =
Conjugate[phi[x1]] phi[x1 + x]/
Sqrt[Abs[phi[x1]]^2 Abs[phi[x1 + x]]^2];
gbar1[x_?NumericQ] :=
1/20 NIntegrate[Evaluate[int[x, x1]], {x1, -10, 10},
Method -> "MonteCarlo", AccuracyGoal -> 2, PrecisionGoal -> 2];
lst1 = Table[{x, Re[gbar1[x]] // Quiet}, {x, -10, 10, 1}];
cor = Interpolation[lst1];

Visualization

Plot[.5 (cor[x] + cor[-x]), {x, -10, 10}]

Note, it is not zero at $$x=\pm 10$$, probably we need accumulate more data.

• Thanks, @AlexTrounev. I will try it in 2D to see If I can get results like arxiv.org/pdf/2107.06680.pdf posted above fig (c). Nov 18, 2022 at 7:17
• Ok, I have to increase the grid points to see the variations. Nov 20, 2022 at 12:32
• @ArghaDebnath Use xmax=60 and plot in the range {x,-20, 20}. Code has been updated. Nov 21, 2022 at 3:52
• Thanks, @AlexTrounev. Nov 21, 2022 at 6:45