Poincaré section of the particle motion near the black hole horizon

Question

I am trying to recreate the Poincaré sections given in Figure 2 of this paper.

For this purpose, I use the following two methods:

Method1: Using ResourceFunction["ClickPoincarePlot2D"]

    f = 2  k  x[t];
lagrangian = -Sqrt[f - x'[t]^2/f - y'[t]^2] - [Omega]^2/
    2  ((x[t] - xc)^2 + y[t]^2);
px1 = D[lagrangian, x'[t]]
py1 = D[lagrangian, y'[t]]
testham = (px1  x'[t] + py1  y'[t]) - lagrangian
k = 1/2; [Omega] = 10; xc = 1;
hamiltonian[x_, px_, y_, py_] = 
  testham /. x[t] -> x /. x'[t] -> px /. y[t] -> y /. y'[t] -> py;
H := Function[{S}, (2  k  S[[1]])/Sqrt[
   2  k  S[[1]] - S[[2]]^2/(2  k  S[[1]]) - S[[4]]^2] + 
   1/2  (S[[1]]^2 + S[[3]]^2 - 2  S[[1]]  xc + xc^2)  [Omega]^2]
{a1, b1, c1, d1} = {D[hamiltonian[x[t], px[t], y[t], py[t]], px[t]], 
  D[hamiltonian[x[t], px[t], y[t], py[t]], 
   py[t]], -D[hamiltonian[x[t], px[t], y[t], py[t]], x[t]], -D[
    hamiltonian[x[t], px[t], y[t], py[t]], y[t]]}; abcd = {x'[t] == 
   a1, px'[t] == c1, y'[t] == b1, py'[t] == d1};

cross = y[t]; recover = {x[t], px[t]};

ResourceFunction[
  "ClickPoincarePlot2D"][abcd, H, 15, t, 10000, cross, recover, \
{PlotStyle -> AbsolutePointSize[1], AspectRatio -> 1, 
  PlotHighlighting -> None, PlotRange -> {{0, 2}, {-20, 20}}}]

Method2: Using algorithm of this answer by @Alex Trounev

psect[{x0_, px0_, y0_, py0_}] := 
 Reap[NDSolve[{abcd, x[0] == x0 + 10^-10, y[0] == y0, px[0] == px0, 
     py[0] == py0, WhenEvent[x[t] == 0, Sow[{y[t], py[t]}]]}, {}, {t, 
     0, 1000}, MaxSteps -> \[Infinity]]][[-1, 1]]
inivalues[E0_, n_, pxmin_, pxmax_] := 
  Module[{x, px, X, y, P, py, Y}, px = RandomReal[{pxmin, pxmax}, n];
   X = Table[0, {n}]; Y = Table[0, {n}]; P = Table[0, {n}];
   Do[{X[[i]], P[[i]], Y[[i]]} = {x, py, y} /. 
       NMinimize[(hamiltonian[x, px[[i]], y, py] - E0)^2, {x, py, 
          y}][[2]];, {i, n}]; Transpose[{X, px, Y, P}]];
iniv = inivalues[15, 100, -0.1, 0.1];
H0[{x_, px_, y_, py_}] = -Sqrt[f - px^2/f - py^2] - [Omega]^2/
    2 ((x - xc)^2 + y^2) + px^2 + py^2;
abcdata = Map[psect, iniv];
ListPlot[abcdata, ImageSize -> Large, Background -> White, 
 Frame -> True]

I tried to alter the cross and recover variables in both methods but to no avail! The first method gives me an empty plot while the second one provides me with errors!

Any help in this regard would be truly beneficial!

Answer 1

There are several problems with your code:

• Hamiltonian construction is not correct (you need to solve for velocities in terms of momenta, not just substitute $x' \to p_x$ and $y' \to p_y$)
• In the linked paper the section is in $(x, p_x)$, while your event is WhenEvent[x[t] == 0, Sow[{y[t], py[t]}]] which is wrong

Note, while solving for $p_x$ and $p_y$ I just took one of the branches. Similarly solving for $p_y$ for given $(x, p_x)$ and $y = 0$, I just assume the last solution is positive. Sections looks close, but you need to check my assumptions.

k = 1/2 ;
w = 10 ; 
c = 1 ;
lagrangian = -Sqrt[2 k x[t] - x'[t]^2/(2 k x[t]) - y'[t]^2] - w^2/2 ((x[t] - c)^2 + y[t]^2) ;
Px = D[lagrangian, x'[t]] ;
Py = D[lagrangian, y'[t]] ;
hamiltonian = ((Px x'[t] + Py y'[t]) - lagrangian) /. (Solve[{px[t] == Px, py[t] == Py}, {x'[t], y'[t]}] // Last) // FullSimplify ;
equations = {
    x'[t] == +D[hamiltonian, px[t]],
    y'[t] == +D[hamiltonian, py[t]],
    px'[t] == -D[hamiltonian, x[t]], 
    py'[t] == -D[hamiltonian, y[t]]
} // Simplify ;
ClearAll[solution] ;
solution[X_, PX_] := Reap@NDSolve[
        Flatten[
            {
                equations, 
                {x[0] == X, px[0] == PX, y[0] == 0, py[0] == initial[X, PX]},
                WhenEvent[y[t] == 0 && y'[t] > 0, Sow[{x[t], px[t]}]]
            }
        ],
        {},
        {t, 0.0, T},
        MaxSteps -> Infinity
] ;
T = 2000.0 ;
energy = 15.0 ;
constraint = (py /. Solve[energy == hamiltonian /. {x[t] -> x, px[t] -> px, y[t] -> 0, py[t] -> py}, py]) // Last ;
ClearAll[initial] ;
initial[x_?NumericQ, px_?NumericQ] := Evaluate[constraint]
ListPlot[Table[solution[X, 0.0] // Last // First, {X, 0.5, 1.5, 0.01}], PlotRange->All, PlotTheme -> "Detailed", PlotStyle -> Directive[{PointSize[Tiny], Black}], AspectRatio -> 1/GoldenRatio, ImageSize->Large, GridLines->None]
energy = 40.0 ;
constraint = (py /. Solve[energy == hamiltonian /. {x[t] -> x, px[t] -> px, y[t] -> 0, py[t] -> py}, py]) // Last ;
ClearAll[initial] ;
initial[x_?NumericQ, px_?NumericQ] := Evaluate[constraint]
ListPlot[Table[solution[X, 0.0] // Last // First, {X, 0.2, 1.5, 0.01}], PlotRange->All, PlotTheme -> "Detailed", PlotStyle -> Directive[{PointSize[Tiny], Black}], AspectRatio -> 1/GoldenRatio, ImageSize->Large, GridLines->None]
energy = 49.5 ;
constraint = (py /. Solve[energy == hamiltonian /. {x[t] -> x, px[t] -> px, y[t] -> 0, py[t] -> py}, py]) // Last ;
ClearAll[initial] ;
initial[x_?NumericQ, px_?NumericQ] := Evaluate[constraint]
ListPlot[Table[solution[X, 0.0] // Last // First, {X, 0.2, 1.5, 0.01}], PlotRange->All, PlotTheme -> "Detailed", PlotStyle -> Directive[{PointSize[Tiny], Black}], AspectRatio -> 1/GoldenRatio, ImageSize->Large, GridLines->None]

Answer 2

To detect more islands, it is better to use the interactive method. Here, I provide an example with the result for the case in which the energy is equal to 40.

PoincareSections[{x0_, px0_, y0_, py0_}, e_, tmax_] := 
 Block[{pdy0, x, px, y, py},
  pdy0 = 
   Sqrt[2500 + e^2 - 100*e*(-1 + x0)^2 + 
      x0*(-10001 - px0^2*x0 + 2500*x0*(6 + (-4 + x0)*x0))]/Sqrt[x0];
  If[Head[N[pdy0]] === Complex, Nothing,
   If[# == {}, {}, First[#]] &@Last[Reap[NDSolve[{
        x'[t] == px[t] x[t] Sqrt[x[t]/(1 + py[t]^2 + px[t]^2 x[t])],
        px'[t] == 
         100 - 100 x[t] - 
          px[t]^2 Sqrt[x[t]/(
           1 + py[t]^2 + px[t]^2 x[t])] - ((1 + py[t]^2) Sqrt[x[t]/(
           1 + py[t]^2 + px[t]^2 x[t])])/(2 x[t]),
        y'[t] == py[t] Sqrt[x[t]/(1 + py[t]^2 + px[t]^2 x[t])],
        py'[t] == -100 y[t],
        x[0] == x0,
        px[0] == px0,
        y[0] == y0,
        py[0] == pdy0,
        WhenEvent[y[t], 
         Sow[{x[t], px[t]}, "LocationMethod" -> "EvenLocator"]]
        },
       {x, px, y, py}, {t, 0, tmax}, MaxSteps -> ∞, 
       MaxStepSize -> 0.1
       ]]]]]

Next, use the code with ClickPane found at: The Poincaré Sections of Yang-Mills-Higgs System:

    style = {PlotStyle -> {{AbsolutePointSize[1], GrayLevel[0], 
             Opacity[0.4]}}, AspectRatio -> 1, ImageSize -> 420};
(*Omit the PlotHighlighting option if you're using a version prior to 13.3.*)
DynamicModule[{f = {}}, 
 Column[{Dynamic[
    MatrixForm@{MousePosition["Graphics", "Mouse not in graphics!"]}],
    Row[{ClickPane[
      Dynamic@ListPlot[f, 
        style], (AppendTo[f, 
         PoincareSections[{#[[1]], #[[2]], 0, 0}, 40, 3000]]) &], 
     Button[Style["Copy orbit list", 12], 
      CopyToClipboard[Iconize[f, "Orbit list"]], 
      ImageSize -> Automatic, Alignment -> Center, 
      Appearance -> {"Palette", "Pressed"}]}, "  "]}, 
  Alignment -> Center]]

The Poincaré section plot looks like this:

style ={PlotStyle -> {{Black, Opacity[0.5], PointSize[0.002]}}, ImageSize -> Large};
ListPlot["Orbit list", style]

Note that you can copy the list of sections by clicking on the button that says Copy orbit list. Afterward, you can save your generated image in PNG format.

Interactively:

The plots of the Poincaré sections for energies 15 and 49.5 respectively look like this:

To export your images with high quality, I recommend the following post: Obtaining better quality in ListPlot.

Answer 3

As mentioned by I.M. the Hamiltonian could be derive with using Legendre transformation - see for example here. Therefore we have

k = 1/2;
w = 10;
c = 1;
lagrangian = -Sqrt[2  k  x[t] - x'[t]^2/(2  k  x[t]) - y'[t]^2] - 
   w^2/2  ((x[t] - c)^2 + y[t]^2);
Px = D[lagrangian, x'[t]];
Py = D[lagrangian, y'[t]];
hamiltonian = ((Px  x'[t] + Py  y'[t]) - 
     lagrangian) /. (Solve[{px[t] == Px, py[t] == Py}, {x'[t], 
       y'[t]}] // First) /. {x[t] -> x, px[t] -> px, y[t] -> y, 
   py[t] -> py}

Last expression we use in two forms

h[x_, px_, y_, py_] := (
   py^2 x)/((1 + py^2 + px^2 x) Sqrt[
    x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
     1 + py^2 + px^2 x)]) + (
   px^2 x^2)/((1 + py^2 + px^2 x) Sqrt[
    x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
     1 + py^2 + px^2 x)]) + Sqrt[
   x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
    1 + py^2 + px^2 x)] + 50  ((-1 + x)^2 + y^2);
H0[{x_, px_, y_, py_}] := (
   py^2 x)/((1 + py^2 + px^2 x) Sqrt[
    x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
     1 + py^2 + px^2 x)]) + (
   px^2 x^2)/((1 + py^2 + px^2 x) Sqrt[
    x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
     1 + py^2 + px^2 x)]) + Sqrt[
   x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
    1 + py^2 + px^2 x)] + 50  ((-1 + x)^2 + y^2);

In this case we can generate initial values using NMiimize as follows

inivalues[E0_, n_, xmin_, xmax_] := 
  Module[{X, x, px, PX, y, P, py, Y}, 
   x = RandomReal[{xmin, xmax}, n];
   PX = Table[0, {n}]; Y = Table[0, {n}]; P = Table[0, {n}];
   Do[{PX[[i]], P[[i]], Y[[i]]} = {px, py, y} /. 
       NMinimize[(h[x[[i]], px, y, py] - E0)^2, {px, py, 
          y}][[2]];, {i, n}]; Transpose[{x, PX, Y, P}]];

For E0=15 we have iniv = inivalues[15, 100, .5, 1.5];

Map[H0, iniv]

(*Out[]= {15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15.}*)

Therefore the energy is saved for every run. Last portion of code is same as in my answer here, also we use {x[t],px[t]} for visualization and y[t] == 0 && y'[t] > 0 for events detection as in a paper linked

{a1, b1, c1, d1} = {D[h[x[t], px[t], y[t], py[t]], px[t]], 
  D[h[x[t], px[t], y[t], py[t]], 
   py[t]], -D[h[x[t], px[t], y[t], py[t]], x[t]], -D[
    h[x[t], px[t], y[t], py[t]], y[t]]}; abcd = {x'[t] == a1, 
  px'[t] == c1, y'[t] == b1, py'[t] == d1};
psect[{x0_, px0_, y0_, py0_}] := 
 Reap[NDSolve[{abcd, x[0] == x0, y[0] == y0, px[0] == px0, 
     py[0] == py0, 
     WhenEvent[y[t] == 0 && y'[t] > 0, Sow[{x[t], px[t]}]]}, {}, {t, 
     0, 1000}, MaxSteps -> 10^6]][[-1, 1]]

Finally we have

abcdata = Map[psect, iniv];
ListPlot[abcdata, ImageSize -> Large, Background -> Black, 
 Frame -> True]

For E0=40 we have picture

And for E0=49.6

Answer 4

The resource function ["ClickPoincarePlot2D"] is also working now, thanks to the new update.

I have also used the ScalarPureFunction by @E. Chan-López presented here to simplify my Hamiltonian into a scalar function form.

ScalarPureFunction[f_, 
    argvars_?
     VectorQ] /; (SameQ[Length[Variables[Level[f, {-1}]]], 
      Length[Flatten[argvars]]] && Length[argvars] > 1) := 
  Module[{positions, heldPart, replacements}, 
   positions = MapIndexed[List, argvars, {-1}];
   heldPart = 
    MapThread[
     Function[
      HoldForm[Part][\[FormalP], 
       Apply[Sequence, Part[SlotSequence@1, 2, All]]]], {positions}];
   replacements = Thread[argvars -> heldPart];
   ReleaseHold[
    Fold[Function[#, #2] &, {\[FormalP]}, {f /. replacements}]]];
ScalarPureFunction[f_][argvars_?VectorQ] := 
  ScalarPureFunction[f, argvars];
k = 1/2;
w = 10;
c = 1;
lagrangian = -Sqrt[2  k  x[t] - x'[t]^2/(2  k  x[t]) - y'[t]^2] - 
   w^2/2  ((x[t] - c)^2 + y[t]^2);
Px = D[lagrangian, x'[t]];
Py = D[lagrangian, y'[t]];
hamiltonian = ((Px  x'[t] + Py  y'[t]) - 
      lagrangian) /. (Solve[{px[t] == Px, py[t] == Py}, {x'[t], 
        y'[t]}] // Last) // FullSimplify;
H := ScalarPureFunction[
  hamiltonian /. {x[t] -> x, px[t] -> px, y[t] -> y, py[t] -> py}, {x,
    px, y, py}]
{a1, b1, c1, d1} = {D[hamiltonian, px[t]], 
  D[hamiltonian, 
   py[t]], -D[hamiltonian, x[t]], -D[hamiltonian, 
    y[t]]}; hameqns = {x'[t] == a1, px'[t] == c1, y'[t] == b1, 
  py'[t] == d1};
cross = y[t]; recover = {x[t], px[t]};

E=15

ResourceFunction[
          "ClickPoincarePlot2D"][hameqns, H, 15, t, 6000, cross, recover,
        {PlotStyle -> AbsolutePointSize[1], AspectRatio -> 1, 
          PlotRange -> {{0, 2}, Automatic}, PlotHighlighting -> None}]

E=40

ResourceFunction
  "ClickPoincarePlot2D"][hameqns, H, 40, t, 6000, cross, recover,
{PlotStyle -> AbsolutePointSize[1], AspectRatio -> 1, 
  PlotRange -> {{0, 2}, Automatic}, PlotHighlighting -> None}]

E=49.5

ResourceFunction
      "ClickPoincarePlot2D"][hameqns, H, 49.5, t, 6000, cross, recover,
    {PlotStyle -> AbsolutePointSize[1], AspectRatio -> 1, 
      PlotRange -> {{0, 2}, Automatic}, PlotHighlighting -> None}]