# ItoProcess for 3 coupled SDEs sourced by 5 Wiener processes

Context

I am trying to solve a stochastic equation corresponding to Vector resonant relaxation (the way orbital planes of stars diffuse near the galactic center, see below).

## Setup

Within the context of some approximation, the orientation vector of the plane of the orbit $$\hat{\mathbf{L}}$$ obeys where $$\mathbf{M}$$ is a stochastic matrix of the form which I have implemented in Mathematica as (assuming $$\Gamma_t=1$$)

Format[em2] = Subscript[e, -2]; Format[em1] = Subscript[e, -1];
Format[e2] = Subscript[e, 2];
Format[e1] = Subscript[e, 1];
Format[e0] = Subscript[e, 0]; Format[u1] = Subscript[u, 1];
Format[u2] = Subscript[u, 2]; Format[u3] = Subscript[u, 3];

M = {{DifferentialD[e2][t],
DifferentialD[em2][t], -DifferentialD[e1][t]}, {DifferentialD[em2][
t], -DifferentialD[e2][t], -DifferentialD[em1][
t]}, {-DifferentialD[e1][t], -DifferentialD[em1][
t], Sqrt DifferentialD[e0][t]}}


Then I can define $$\hat{ \mathbf{L}}$$ as

L = #[t] & /@ {u1, u2, u3};
eqn = \[DifferentialD]#[t] & /@ {u1, u2, u3} == Cross[L, M.L] // Apart // Thread;
eqn // TableForm So that the ItoProcess obeys

proc = ItoProcess[
eqn, {u1[t], u2[t], u3[t]}, {{u1, u2, u3}, {1, 1/2, 1/3}}, t,
{e2 \[Distributed] WienerProcess[],
em2 \[Distributed] WienerProcess[],
e1 \[Distributed] WienerProcess[],
em1 \[Distributed] WienerProcess[],
e0 \[Distributed] WienerProcess[]}]


I can then compute a set of paths as

path = RandomFunction[proc, {0., 1., 0.01}, Method -> "StochasticRungeKutta"]
ListLinePlot[path] Question

Why is it not returning random paths?

Let me apologize in advance for what is likely to be a stupid mistake on my part.

Astrophysical context

In the vicinity of a supermassive black hole, stars move on nearly Keplerian orbits. Yet, because of the enclosed stellar mass and relativistic corrections, the potential slightly deviates from the Keplerian one, which causes the stellar orbits to precess. Similarly, as a result of the finite number of stars, the mutual gravitational torques between pairs of stars also drive a rapid reshuffling of the stars' orbital orientations. Overall, the combination of these two effects leads to a stochastic evolution of stellar orbital angular momentum vectors, through a process named 'resonant relaxation'.

The typographic error is in the definition of the matrix

M = {{DifferentialD[e2[t]],
DifferentialD[em2[t]], -DifferentialD[e1[t]]}, {DifferentialD[
em2[t]],
-DifferentialD[e2[t]], -DifferentialD[em1[t]]}, {-DifferentialD[
e1[t]], -DifferentialD[em1[t]], Sqrt DifferentialD[e0[t]]}}


i.e. it should involve DifferentialD[e2[t]] or \[DifferentialD]e2[t] not DifferentialD[e2][t]. Since for this problem the noise needs to be correlated (see this question), it is in fact simpler to rely on NDSolve together with GaussianRandomField

nn = 1024*2; tmax = 2*25; Amp = 5;
noise = Interpolation[#, Method -> "Spline"][t] & /@
Table[Transpose@{tmax Range[0, nn - 1]/(nn - 1) // N,
GaussianRandomField[nn, 1, Function[k, k^-0.5]] Amp}, {5}];
eqn2 = eqn /. {\[DifferentialD]e0[t] ->
noise[], \[DifferentialD]e1[t] ->
noise[], \[DifferentialD]em1[t] ->
noise[], \[DifferentialD]e2[t] ->
noise[], \[DifferentialD]em2[t] ->
noise[]} /.
{\[DifferentialD]u1[t] ->
u1'[t], \[DifferentialD]u2[t] ->
u2'[t], \[DifferentialD]u3[t] -> u3'[t]};
eqn2 = Join[eqn2, {u1 == 1, u2 == 0, u3 == 0}];
sol[t_] = NDSolveValue[eqn2, {u1[t], u2[t], u3[t]}, {t, 0, tmax}]


Then

Show[ParametricPlot3D[ sol[t], {t, 0, tmax}, BoxRatios -> {1, 1, 1}],
Graphics3D[{Opacity[0.2], Sphere[]}], Boxed -> False, Axes -> False] One can represent explicitly the orbital planes via this function

circle3D[centre_: {0, 0, 0}, radius_: 1, normal_: {0, 0, 1},
angle_: {0, 2 Pi}] :=
Composition[Tube,
Map[RotationTransform[{{0, 0, 1}, normal}, centre], #] &,
Map[Append[#, Last@centre] &, #] &,
Append[DeleteDuplicates[Most@#], Last@#] &, Level[#, {-2}] &,
MeshPrimitives[#, 1] &, DiscretizeRegion, If][
First@Differences@angle >= 2 Pi, Circle[Most@centre, radius],

Show[Graphics3D[{Thick, 