ClearAll["Global`*"]
p1 = {l Sin[\[Theta][t]], -l Cos[\[Theta][t]]};
p2 = p1 + {l Sin[\[Phi][t]], -l Cos[\[Phi][t]]};
T = FullSimplify[1/2 m Norm[D[p1, t]]^2 + 1/2 m Norm[D[p2, t]]^2, Assumptions -> {m > 0, l > 0, \[Theta][t] > 0, \[Theta]'[t] > 0, \[Phi]'[t] > 0, \[Phi][t] > 0}];
V = FullSimplify[m g p1[[2]] + m g p2[[2]], Assumptions -> {m > 0 , l > 0}];
L = FullSimplify[T - V];
cmomentum\[Theta][t_] := FullSimplify[D[L, \[Theta]'[t]] ];
cmomentum\[Phi][t_] := FullSimplify[D[L, \[Phi]'[t]] ];
momentumsub = Solve[{p\[Theta][t] == cmomentum\[Theta][t], p\[Phi][t] == cmomentum\[Phi][t]}, {\[Theta]'[t], \[Phi]'[t]}];
H = FullSimplify[(p\[Theta][t]*\[Theta]'[t] + p\[Phi][t]*\[Phi]'[t] - L) /. momentumsub][[1]];

tmax = 100;
m = 1; g = 1; l = 1;
\[Theta]0 = \[Pi]/2; p\[Theta]0 = 0; \[Phi]0 = \[Pi]/2; p\[Phi]0 = 0;
n = 5000;

eqs = {\[Theta]'[t] == D[H, p\[Theta][t]], p\[Theta]'[t] == -D[H, \[Theta][t]], \[Phi]'[t] == D[H, p\[Phi][t]], p\[Phi]'[t] == -D[H, \[Phi][t]]};

initial = RandomVariate[MultinormalDistribution[{\[Theta]0, p\[Theta]0, \[Phi]0, p\[Phi]0}, 0.0001 IdentityMatrix[4]], n];
sols = Table[
   NDSolve[{eqs, \[Theta][0] == initial[[j, 1]], p\[Theta][0] == initial[[j, 2]], \[Phi][0] == initial[[j, 3]], p\[Phi][0] == initial[[j, 4]]}, {p\[Theta][t], \[Theta][t], p\[Phi][t], \[Phi][t]}, {t, 0, tmax}][[1]], {j, 1, n}];
npend = 100;
frames = Table[
   Grid[{{
      Graphics[{
        Thick, Disk[{0, 0}, 0.05],
        Opacity[0.5],
        Table[ Disk[p1 /. sols[[j]], 0.1], {j, 1, npend}], 
        Table[Disk[p2 /. sols[[j]], 0.1], {j, 1, npend}],
        Table[ Line[{{0, 0} + 0.05 Normalize[p1], (1 - 0.1) p1} /. sols[[j]]], {j, 1, npend}],
        Table[ Line[{p1 + 0.1 Normalize[p2 - p1], p2 - 0.1 Normalize[p2 - p1]} /. sols[[j]]], {j, 1, npend}]
        }, PlotRange -> 2.2 {{-l, l}, {-l, l}}]
      ,
      Grid[{{
         Graphics[{Opacity[0.2],
           Black, PointSize[0.01], 
           Point[{Mod[\[Theta][t], 2 \[Pi], -\[Pi]], p\[Theta][t]} /. sols]
           }, PlotRange -> (1.2*{{-\[Pi], \[Pi]}, {-\[Pi], \[Pi]}}), AspectRatio -> 1, Axes -> False, Frame -> True, FrameTicks -> None, FrameLabel -> {"\[Theta]", "\!\(\*SubscriptBox[\(p\), \(\[Theta]\)]\)"}, LabelStyle -> {Black, Bold, FontSize -> 14}],
         Graphics[{Opacity[0.2],
           Black, PointSize[0.01], 
           Point[{Mod[\[Phi][t], 2 \[Pi], -\[Pi]], p\[Theta][t]} /. sols]
           }, PlotRange -> (1.2*{{-\[Pi], \[Pi]}, {-\[Pi], \[Pi]}}), AspectRatio -> 1, Axes -> False, Frame -> True, FrameTicks -> None, FrameLabel -> {"\[Phi]", "\!\(\*SubscriptBox[\(p\), \(\[Theta]\)]\)"},           LabelStyle -> {Black, Bold, FontSize -> 14}]
         }, {
         Graphics[{Opacity[0.2],
           Black, PointSize[0.01], 
           Point[{Mod[\[Theta][t], 2 \[Pi], -\[Pi]], p\[Phi][t]} /. sols]
           }, PlotRange -> (1.2*{{-\[Pi], \[Pi]}, {-\[Pi], \[Pi]}}), AspectRatio -> 1, Axes -> False, Frame -> True, FrameTicks -> None, FrameLabel -> {"\[Theta]", "\!\(\*SubscriptBox[\(p\), \(\[Phi]\)]\)"},           LabelStyle -> {Black, Bold, FontSize -> 14}],
         Graphics[{Opacity[0.2],
           Black, PointSize[0.01], 
           Point[{Mod[\[Phi][t], 2 \[Pi], -\[Pi]], p\[Phi][t]} /. sols]
           }, PlotRange -> (1.2*{{-\[Pi], \[Pi]}, {-\[Pi], \[Pi]}}), AspectRatio -> 1, Axes -> False, Frame -> True, FrameTicks -> None, FrameLabel -> {"\[Phi]", "\!\(\*SubscriptBox[\(p\), \(\[Phi]\)]\)"},           LabelStyle -> {Black, Bold, FontSize -> 14}]
         }}]
      }}]
   , {t, 0, 45, 0.25}];