In this simple computational essay, I investigate the dynamics of a Bose-Einstein condensate (BEC) governed by the Gross-Pitaevskii equation under the influence of a time-dependent perturbation. The perturbation consists of an attractive localized potential traversing the condensate along various trajectories. Numerical simulations of the dynamical equation reveal the emergence of some sort of wake patterns in response to the perturbation. While this is a simple study, I believe it hints at deeper concepts and applications yet to be explored.
RELATED READING: G W Stagg et al 2014 J. Phys. B: At. Mol. Opt. Phys. 47 095304
Berry, M V, 2001 'Why are special functions special?', Physics Today, April, pp11-12
RELATED READING: G W Stagg et al 2014 J. Phys. B: At. Mol. Opt. Phys. 47 095304
Berry, M V, 2001 'Why are special functions special?', Physics Today, April, pp11-12
This exploration was inspired by discussions with Sir Michael Berry during the Wolfram Innovator Award 2025 meetings. Also, Michael Trott suggested to examine trajectories beyond straight lines.
Line
Line
Trajectory of perturbation:
In[]:=
With{initial=DensityPlot[Exp[-(x^2+y^2)],{x,-15,15},{y,-10,10},PlotRange->All,PlotPoints->50,ColorFunction->"SunsetColors"]},ManipulateShowinitial,Quiet@ContourPlotExp[-(x-7τ+7)^2-],{x,-15,15},{y,-10,10},PlotRange->All,ContourShading->None,PlotPoints->50,ContourStyle->Red,ParametricPlot[{7u-7,0},{u,0,3},PlotStyle->Directive[Dotted,Green]],PlotLabel->"Perturbation wrt the initial wavefunction",{τ,0,3},SaveDefinitions->True
2
(y)
Hamiltonian:
In[]:=
ℋ=SchrodingerPDEComponent{Ψ[t,x,y],t,{x,y}},<|"ReducedPlanckConstant"->1,"SchrodingerPotential"->3Ψ[t,x,y]Conjugate[Ψ[t,x,y]]+.5(x^2+y^2)-29Exp[-(x-7t+7)^2-]|>;
2
(y)
Initial state:
In[]:=
ic=Ψ[0,x,y]==Exp[-(x^2+y^2)];
Solve dynamical equation numerically:
In[]:=
sol=NDSolveValue[{ℋ==0,ic},Ψ,{t,0,3},{x,-20,20},{y,-10,10},Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->"MaxCellMeasure"->.2}}}];//AbsoluteTiming
Visualize results:
In[]:=
AnimationVideo[Show[DensityPlot[Abs[sol[t,x,y]],{x,-20,20},{y,-10,10},PlotPoints150,ColorFunction(ColorData["SunsetColors"][Sqrt[#]]&),PlotRange->All],Graphics[{Green,PointSize[Medium],Point[{7t-7,0}]}]],{t,0,3}]
Get a snapshots of evolution:
In[]:=
imgs=Show[DensityPlot[Abs[sol[#,x,y]],{x,-20,20},{y,-10,10},PlotPoints150,ColorFunction(ColorData["SunsetColors"][Sqrt[#]]&),PlotRange->All],Graphics[{Red,PointSize[Medium],Point[{7#-7,0}]}]]&/@Range[0,3,.05];
Export as GIF:
In[]:=
Export[NotebookDirectory[]<>"Line.gif",imgs,"DisplayDurations"->.3,AnimationRepetitions->100]
Circle
Circle
Radius:
In[]:=
R=1;
Trajectory of perturbation:
With[{initial=DensityPlot[Exp[-(x^2+y^2)],{x,-3,3},{y,-3,3},PlotRange->All,PlotPoints->50,ColorFunction->"SunsetColors"]},Manipulate[Show[initial,Quiet@ContourPlot[Exp[-(x-RCos[2πτ])^2-(y-RSin[2πτ])^2],{x,-3,3},{y,-3,3},PlotRange->All,ContourShading->None,PlotPoints->100,ContourStyle->Red],ParametricPlot[R{Cos[2πu],Sin[2πu]},{u,0,1},PlotStyle->Directive[Dotted,Green]],PlotLabel->"Perturbation wrt the initial wavefunction"],{τ,0,3},SaveDefinitions->True]]
Hamiltonian:
In[]:=
ℋ=SchrodingerPDEComponent[{Ψ[t,x,y],t,{x,y}},<|"ReducedPlanckConstant"->1,"SchrodingerPotential"->3Ψ[t,x,y]Conjugate[Ψ[t,x,y]]+.5(x^2+y^2)-29Exp[-(x-RCos[2πt])^2-(y-RSin[2πt])^2]|>];
Initial state:
In[]:=
ic=Ψ[0,x,y]==Exp[-(x^2+y^2)];
Solve the dynamical equation numerically:
In[]:=
sol=NDSolveValue[{ℋ==0,ic},Ψ,{t,0,2},{x,-15,15},{y,-15,15},Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->"MaxCellMeasure"->.1}}}];//AbsoluteTiming
Visualize results:
In[]:=
AnimationVideo[Show[DensityPlot[Abs[sol[t,x,y]],{x,-10,10},{y,-10,10},PlotPoints150,ColorFunction(ColorData["SunsetColors"][Sqrt[#]]&),PlotRange->All],Graphics[{Green,Point[R{Cos[2πt],Sin[2πt]}]}]],{t,0,2}]
A series of snapshots:
In[]:=
imgs=Table[Show[DensityPlot[Abs[sol[t,x,y]],{x,-10,10},{y,-10,10},PlotPoints150,ColorFunction(ColorData["SunsetColors"][Sqrt[#]]&),PlotRange->All],Graphics[{Green,Point[R{Cos[2πt],Sin[2πt]}]}]],{t,0,2,.05}];
Export as GIF:
In[]:=
Export[NotebookDirectory[]<>"Circle.gif",imgs,"DisplayDurations"->.3,AnimationRepetitions->100]
Spiral
Spiral
Trajectory of perturbation
In[]:=
With[{initial=DensityPlot[Exp[-(x^2+y^2)],{x,-3,3},{y,-3,3},PlotRange->All,PlotPoints->50,ColorFunction->"SunsetColors"]},Manipulate[Show[initial,Quiet@ContourPlot[Exp[-(x-τCos[2πτ])^2-(y-τSin[2πτ])^2],{x,-3,3},{y,-3,3},PlotRange->All,ContourShading->None,PlotPoints->100,ContourStyle->Red],ParametricPlot[{uCos[2πu],uSin[2πu]},{u,0,2},PlotStyle->Directive[Dotted,Green]],PlotLabel->"Perturbation wrt the initial wavefunction"],{τ,0,2},SaveDefinitions->True]]
Hamiltonian:
In[]:=
ℋ=SchrodingerPDEComponent[{Ψ[t,x,y],t,{x,y}},<|"ReducedPlanckConstant"->1,"SchrodingerPotential"->3Ψ[t,x,y]Conjugate[Ψ[t,x,y]]+.5(x^2+y^2)-29Exp[-(x-tCos[2πt])^2-(y-tSin[2πt])^2]|>];
Initial state:
In[]:=
ic=Ψ[0,x,y]==Exp[-(x^2+y^2)];
Solve the dynamical equation numerically:
In[]:=
sol=NDSolveValue[{ℋ==0,ic},Ψ,{t,0,2},{x,-15,15},{y,-15,15},Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->"MaxCellMeasure"->.1}}}];//AbsoluteTiming
Visualize results:
In[]:=
AnimationVideo[Show[DensityPlot[Abs[sol[t,x,y]],{x,-10,10},{y,-10,10},PlotPoints150,ColorFunction(ColorData["SunsetColors"][Sqrt[#]]&),PlotRange->All],Graphics[{Green,Point[{tCos[2πt],tSin[2πt]}]}]],{t,0,2}]
A series of snapshots:
In[]:=
imgs=Table[Show[DensityPlot[Abs[sol[t,x,y]],{x,-10,10},{y,-10,10},PlotPoints150,ColorFunction(ColorData["SunsetColors"][Sqrt[#]]&),PlotRange->All],Graphics[{Green,Point[{tCos[2πt],tSin[2πt]}]}]],{t,0,2,.05}];
Export as GIF:
In[]:=
Export[NotebookDirectory[]<>"Spiral.gif",imgs,"DisplayDurations"->.3,AnimationRepetitions->100]
Bowtie
Bowtie
Trajectory of perturbation
In[]:=
With[{initial=DensityPlot[Exp[-(x^2+y^2)],{x,-3,3},{y,-3,3},PlotRange->All,PlotPoints->50,ColorFunction->"SunsetColors"]},Manipulate[Show[initial,Quiet@ContourPlot[Exp[-(x-1.5Sin[2πτ])^2-(y-1.5Sin[4πτ])^2],{x,-3,3},{y,-3,3},PlotRange->All,ContourShading->None,PlotPoints->100,ContourStyle->Red],ParametricPlot[1.5{Sin[2πu],Sin[4πu]},{u,0,2},PlotStyle->Directive[Dotted,Green]],PlotLabel->"Perturbation wrt the initial wavefunction"],{τ,0,2},SaveDefinitions->True]]
Hamiltonian:
In[]:=
ℋ=SchrodingerPDEComponent[{Ψ[t,x,y],t,{x,y}},<|"ReducedPlanckConstant"->1,"SchrodingerPotential"->3Ψ[t,x,y]Conjugate[Ψ[t,x,y]]+.5(x^2+y^2)-29Exp[-(x-1.5Sin[2πt])^2-(y-1.5Sin[4πt])^2]|>];
Initial state:
In[]:=
ic=Ψ[0,x,y]==Exp[-(x^2+y^2)];
Solve the dynamical equation numerically:
In[]:=
sol=NDSolveValue[{ℋ==0,ic},Ψ,{t,0,2},{x,-15,15},{y,-15,15},Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->"MaxCellMeasure"->.1}}}];//AbsoluteTiming
Visualize results:
In[]:=
AnimationVideo[Show[DensityPlot[Abs[sol[t,x,y]],{x,-10,10},{y,-10,10},PlotPoints150,ColorFunction(ColorData["SunsetColors"][Sqrt[#]]&),PlotRange->All],Graphics[{Green,Point[{Sin[2πt],Sin[2π2t]}]}]],{t,0,2}]
A series of snapshots:
In[]:=
imgs=Table[Show[DensityPlot[Abs[sol[t,x,y]],{x,-10,10},{y,-10,10},PlotPoints150,ColorFunction(ColorData["SunsetColors"][Sqrt[#]]&),PlotRange->All],Graphics[{Green,Point[{Sin[2πt],Sin[2π2t]}]}]],{t,0,2,.05}];
Export as GIF:
In[]:=
Export[NotebookDirectory[]<>"BowTie.gif",imgs,"DisplayDurations"->.3,AnimationRepetitions->100]
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Quantum Wakes: Bose-Einstein condensate under time-dependent perturbations
by Mohammad Bahrami
Wolfram Community, STAFF PICKS, December 2, 2024
https://community.wolfram.com/groups/-/m/t/3331281
by Mohammad Bahrami
Wolfram Community, STAFF PICKS, December 2, 2024
https://community.wolfram.com/groups/-/m/t/3331281