12
$\begingroup$

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.

enter image description here

 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}}]

enter image description here

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

enter image description here

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.

enter image description here

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]

enter image description here

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]

enter image description here

$\endgroup$
11
  • $\begingroup$ What's wrong with the error message you get? That needs to be fixed first. $\endgroup$
    – user21
    Nov 3, 2022 at 10:28
  • $\begingroup$ 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. $\endgroup$
    – Chris K
    Nov 3, 2022 at 13:28
  • $\begingroup$ Thanks for the reference @Chris K. $\endgroup$ Nov 3, 2022 at 14:05
  • $\begingroup$ @ArghaDebnath I see you've edited your question to use ItoProcess now. Did it work or are you still having problems? $\endgroup$
    – Chris K
    Nov 4, 2022 at 13:02
  • $\begingroup$ No, @Chris K it is running for a long time with no answer. $\endgroup$ Nov 4, 2022 at 16:08

1 Answer 1

4
+50
$\begingroup$

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]

Figure 1

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}]

Figure 2

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

$\endgroup$
4
  • $\begingroup$ 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). $\endgroup$ Nov 18, 2022 at 7:17
  • $\begingroup$ Ok, I have to increase the grid points to see the variations. $\endgroup$ Nov 20, 2022 at 12:32
  • $\begingroup$ @ArghaDebnath Use xmax=60 and plot in the range {x,-20, 20}. Code has been updated. $\endgroup$ Nov 21, 2022 at 3:52
  • $\begingroup$ Thanks, @AlexTrounev. $\endgroup$ Nov 21, 2022 at 6:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.