CITE THIS NOTEBOOK: Dynamics of stochastic-constrained particles by Tao Guo. Wolfram Community FEB 21 2023.
ORIGINAL ARTICLE: Guo, T. Dynamics of stochastic-constrained particles. Scientific Reports, Nature 13, 2759 (2023). https://doi.org/10.1038/s41598-023-29940-y.
Article Abstract: Prior studies have focused on the overall behavior of randomly moving particle swarms. However, the characteristics of the stochastic-constrained particles that form ubiquitously within these swarms remain oblivious. This study demonstrates a generalized diffusion equation for stochastic-constrained particles that considers the velocity and location aggregation effects observed from their parent particle swarm (i.e., a completely random particle swarm). This equation can be approximated as the form of Schrödinger equation in the microcosmic case (low relative density) and describe the dynamics of the total mass distribution in the macrocosmic case (high relative density). The predicted density distribution of the particle swarm in the stable aggregation state is consistent with the total mass distribution of massive, relaxed galaxy clusters (at least in the range of
r<
r
s
), preventing cuspy problems in the empirical Navarro–Frenk–White profile. This study opens a window to observe the dynamics of stochastic-constrained particles from a third perspective, from which the aggregation effect of particles without gravitation can be saw.
NOTE:
1. The "Euclid Math One" regular and bold fonts are needed to display the contents correctly in this Notebook.
2. If there is no special case, the Mathematica code starts with gray "
In[]:=
" and is bold by default according to Mathematica's rules.
Biased stochastic particles can exhibit various forms of bias. If the particles’ motion shows aggregation of velocity direction, a special relativity-like effect arises within the particle swarm. If it shows positional aggregation, an effect equivalent to gravitation is produced. The main text deduces the dynamic equation that describes the motion of such particles:

∂M
∂t
=-
ℏ
2
Γ
2m
M

[
2
∇
M-
2
T
(M)]
, where
2
T
(M)=
2
∂M
∂x
+
2
∂M
∂y
+
2
∂M
∂z
. One of its stable states is
2
∇
M=
2
T
(M)
. This stable state corresponds to a state of large, massive and ancient galaxies or galaxy clusters solely governed by gravitation. Coincidentally, the analytical solution of this stable state equation (in two and three dimensions) can be obtained through Mathematica. It should be noted that the value of this analytical solution can only be calculated in the complex domain. As described in the main text, the norm of this value is proportional to the mass density of the galaxies or galaxy clusters. The results predicted by this equation are consistent with observation on a scale of three orders of magnitude. Part 3 of the Supplementary Information presents the specific process of obtaining the analytical solution for the stable state equation. Part 4 of the Supplementary Information shows the method for generating each figure in the main text, which includes schematic diagrams of the process for deducing the dynamic equation and some figures of the particle motion represented by the equation as it reaches its final state. However, unfortunately, this dynamic equation is difficult to solve, and the main text can only obtain numerical solutions near the initial stage. It is difficult to obtain solutions close to the stable state.

Solving Process of Eq. 38 in the Main Text

To solve the partial differential equation Eq. 38 in the main text, it is assumed that the system is spherically symmetric because it is isotropic at a huge scale. Therefore, we make the conversion from rectangular to spherical coordinates (note that φ is used to denote the azimuthal angle, whereas θ is used to denote the polar angle), namely, x = r sin θ cos φ, y = r sin θ sin φ and z = cos θ.
In the case of spherical symmetry, the change of function M(r) does not depend on θ and φ, but is related to r. Therefore, after the coordinate transformation, and the first and the second derivatives are obtained, to omit the terms that depends on angles θ and φ, we can obtain (subject to the character limitation of Mathematica, ℳ is used instead of M in the code cell; the same is done below):
In[]:=
Simplify
2
r
D[ℳ[r],{r,1}]+D[ℳ[r],{r,2}]-
2
(D[ℳ[r],{r,1}])
(
2
(Sin[θ]Cos[φ])
+
2
(Sin[θ]Sin[φ])
+
2
(Cos[θ])
)
Out[]=
′′
ℳ
(r)-
2
′
ℳ
(r)
+
2
′
ℳ
(r)
r
To solve the abovementioned differential equation under the boundary condition M
(
r
e
)=0
.
In[]:=
DSolve
′′
ℳ
[r]-
2
(
′
ℳ
[r])
+
2
r
′
ℳ
[r]0,ℳ[re]0,ℳ[r],r
Out[]=
{{ℳ(r)log(r)-log(1+
c
1
r)-log(re)+log(1+
c
1
re)}}
Suppose another boundary condition is M
(
r
c
)=
M
c
, then
In[]:=
r=rc;​​Solve[Log[r]-Log[1+c1r]-Log[re]+Log[1+c1re]ℳc,c1]
Out[]=
c1
rc-re
ℳc

rcre(
ℳc

-1)

Therefore, the solution of the above differential equation is as follows:
In[]:=
Clear["Global`*"];​​c1=
rc-re
ℳc

rcre(
ℳc

-1)
;​​Simplify[Log[r]-Log[1+c1r]-Log[re]+Log[1+c1re]]
Out[]=
-log
r(rc-re
ℳc

)
rcre(
ℳc

-1)
+1+log(r)+log
ℳc

(rc-re)
rc(
ℳc

-1)
-log(re)
To restore the above solution in spherical to the solution in 3-dimensional rectangular coordinates, then
In[]:=
r=
2
x
+
2
y
+
2
z
;​​FullSimplify-Log
r(rc-re
ℳc

)
rcre(
ℳc

-1)
+1+Log[r]+Log
ℳc

(rc-re)
rc(
ℳc

-1)
-Log[re],Assumptionsre>rc>0
Out[]=
-log
(rc-re
ℳc

)
2
x
+
2
y
+
2
z
ℳc

-1
+rcre+log
ℳc

(rc-re)
ℳc

-1
+
1
2
log(
2
x
+
2
y
+
2
z
)
To verify the above results:
In[]:=
ℳ[x,y,z]:=-Log
(c-re
ℳc

)
2
x
+
2
y
+
2
z
ℳc

-1
+rcre+Log
ℳc

(rc-re)
ℳc

-1
+
1
2
Log[
2
x
+
2
y
+
2
z
];​​FullSimplify
2
∂
ℳ(x,y,z)
∂
2
x
+
2
∂
ℳ(x,y,z)
∂
2
y
+
2
∂
ℳ(x,y,z)
∂
2
z
-
2
∂ℳ(x,y,z)
∂x
-
2
∂ℳ(x,y,z)
∂y
-
2
∂ℳ(x,y,z)
∂z

Out[]=
0
Therefore, the above equation is the solution of Eq. 38 in the main text (only when
Im[
M
0
]∈[-π,π)
and the principal values of arguments are taken in the calculation process).
Similarly, the 2-dimensional case can also be solved.
In[]:=
Clear["Global`*"];​​SimplifyD[ℳ[r],{r,2}]+
1
r
D[ℳ[r],{r,1}]-
2
(D[ℳ[r],{r,1}])

Out[]=
′′
ℳ
(r)-
2
′
ℳ
(r)
+
′
ℳ
(r)
r
In[]:=
DSolve
′′
ℳ
[r]-
2
′
ℳ
[r]
+
′
ℳ
[r]
r
0,ℳ[re]0,ℳ[r],r
Out[]=
{{ℳ(r)log(-log(re)+
c
1
)-log(-log(r)+
c
1
)}}
In[]:=
r=rc;​​Solve[Log[-Log[re]+c1]-Log[-Log[r]+c1]ℳc,c1]
Out[]=
c1
ℳc

log(rc)-log(re)
ℳc

-1

In[]:=
Clear["Global`*"];​​c1=
ℳc

Log[rc]-Log[re]
ℳc

-1
;​​Simplify[Log[-Log[re]+c1]-Log[-Log[r]+c1]]
Out[]=
log
ℳc

(log(rc)-log(re))
ℳc

-1
-log
ℳc

log(rc)-log(re)
ℳc

-1
-log(r)
In[]:=
r=
2
x
+
2
y
;​​FullSimplifyLog
ℳc

(Log[rc]-Log[re])
ℳc

-1
-Log
ℳc

Log[rc]-Log[re]
ℳc

-1
-Log[r],Assumptionsre>rc>0
Out[]=
log
ℳc

log
rc
re

ℳc

-1
-log
log
rc
re

ℳc

-1
+log(rc)-
1
2
log(
2
x
+
2
y
)
In[]:=
ℳ[x,y]:=Log
ℳc

Log
rc
re

ℳc

-1
-Log
Log
rc
re

ℳc

-1
+Log[rc]-
1
2
Log[
2
x
+
2
y
];​​FullSimplify
2
∂
ℳ(x,y)
∂
2
x
+
2
∂
ℳ(x,y)
∂
2
y
-
2
∂ℳ(x,y)
∂x
-
2
∂ℳ(x,y)
∂y

Out[]=
0
To verify the above conclusion, the results of analytical solution and the numerical solution under the same conditions are plotted (This code takes approximately 32 seconds):
In[]:=
Clear["Global`*"];​​ℳa[x_,y_]:=Log
ℳc

Log
rc
re

ℳc

-1
-Log
Log
rc
re

ℳc

-1
+Log[rc]-
1
2
Log[
2
x
+
2
y
];​​rc=
4
100
;​​re=4;​​ℳc=1+2;​​Ω=ImplicitRegion[
2
rc
≤
2
x
+
2
y
≤
2
re
,{x,y}];​​G1=ShowPlot3DNorm[ℳa[x,y]],{x,y}∈Ω,PlotRange0,
8
,ColorFunction(Hue[0.65,#3]&),MeshStyleNone,BoundaryStyleNone,PlotPoints300,AxesLabelStyle["x ",Italic],Style["y",Italic],RotateStyle["Density "],
π
2
,AxesStyleDirective[Black,FontFamily"Arial",FontSize15],TicksStyleBlack,BoxStyleDirective[Black,Thickness0.0018],BoxRatiosAutomatic,ViewPoint{15,-26,16},EpilogText[Style["a",15,FontFamily"Arial",Bold,Black],{-0.07,0.92},{-1,1}],TableΩ1=ImplicitRegion
9
100
≤
2
x
+
2
i
≤16,{x};If
2
i
≤
9
100
,xx=
9
100
-
2
i
,xx=0;ParametricPlot3D{x,i,Norm[ℳa[x,i]]},{x}∈Ω1,PlotStyleThickness[0.0018],PlotPoints300,ColorFunctionGrayLevel0.4,1-#3×
Norm[ℳa[xx,i]]
Normℳa0,
3
10

&,{i,-3.5,3.5,0.5},TableΩ1=ImplicitRegion
9
100
≤
2
j
+
2
y
≤16,{y};If
2
j
≤
9
100
,yy=
9
100
-
2
j
,yy=0;ParametricPlot3D{j,y,Norm[ℳa[j,y]]},{y}∈Ω1,PlotStyleThickness[0.0018],PlotPoints300,ColorFunctionGrayLevel0.4,1-#3×
Norm[ℳa[j,yy]]
Normℳa0,
3
10

&,{j,-3.5,3.5,0.5},ParametricPlot3D[{4Cos[ϕ],4Sin[ϕ],0},{ϕ,0,2π},PlotStyleDirective[Gray,Thickness[0.0018]],PlotPoints300];​​Needs["NDSolve`FEM`"];​​mesh=ToElementMeshΩ,MeshRefinementFunctionFunction{vertices,area},area>
3
100000
1
10
+80Norm[Mean[vertices]];​​ℳn=NDSolveValue
2
∂
u(x,y)
∂
2
x
+
2
∂
u(x,y)
∂
2
y
-
2
∂u(x,y)
∂x
-
2
∂u(x,y)
∂y
0,DirichletCondition[u[x,y]ℳc,
2
x
+
2
y

2
rc
],DirichletCondition[u[x,y]0,
2
x
+
2
y

2
re
],u,{x,y}∈mesh;​​G2=ShowPlot3DNorm[ℳn[x,y]],{x,y}∈mesh,PlotRange0,
8
,ColorFunction(Hue[0.65,#3]&),MeshStyleNone,BoundaryStyleNone,AxesLabelStyle["x ",Italic],Style["y",Italic],RotateStyle["Density "],
π
2
,AxesStyleDirective[Black,FontFamily"Arial",FontSize15],TicksStyleBlack,BoxStyleDirective[Black,Thickness0.002],BoxRatiosAutomatic,ViewPoint{15,-26,16},EpilogText[Style["b",15,FontFamily"Arial",Bold,Black],{-0.07,0.92},{-1,1}],TableΩ2=ImplicitRegion
9
100
≤
2
x
+
2
i
≤16,{x};If
2
i
≤
9
100
,xx=
9
100
-
2
i
,xx=0;ParametricPlot3D{x,i,Norm[ℳn[x,i]]},{x}∈Ω2,PlotStyleThickness[0.0018],PlotPoints300,ColorFunctionGrayLevel0.4,1-#3×
Norm[ℳn[xx,i]]
Normℳn0,
3
10

&,{i,-3.5,3.5,0.5},TableΩ2=ImplicitRegion
9
100
≤
2
j
+
2
y
≤16,{y};If
2
j
≤
9
100
,yy=
9
100
-
2
j
,yy=0;ParametricPlot3D{j,y,Norm[ℳn[j,y]]},{y}∈Ω2,PlotStyleThickness[0.0018],PlotPoints300,ColorFunctionGrayLevel0.4,1-#3×
Norm[ℳn[j,yy]]
Normℳn0,
3
10

&,{j,-3.5,3.5,0.5},ParametricPlot3D[{4Cos[ϕ],4Sin[ϕ],0},{ϕ,0,2π},PlotStyleDirective[Gray,Thickness[0.0018]],PlotPoints300];​​s1=GraphicsRow[{G1,G2},ImageSize500,SpacingsScaled[-0.06]];​​Pane[s1,{500,200},ImageMargins{{50,-30},{-18,-25}}]

Figures Used in the Main Text