# The Poincaré Sections of Yang-Mills-Higgs System [closed]

In the papers Quantum Chaos in a Yang-Mills-Higgs System and DIFFERENT FACETS OF CHAOS IN QUANTUM MECHANICS, the authors have plotted three Poincaré sections for different values of v.

I have been trying to reproduce these sections; however, they are not coming exactly. The following is my code, which I borrowed from this answer.

Hamiltonian is given by

H[px_, py_, x_, y_] := 1/2 (px^2 + py^2) + g^2*v^2*(x^2 + y^2) + 1/2*g^2*x^2*y^2


Values of parameters

   g = 1; v = 1.0(*1.1*)(*1.2*);


Hamiltonian EOMs

{a1, b1, c1, d1} = {D[H[px[t], py[t], x[t], y[t]], px[t]],
D[H[px[t], py[t], x[t], y[t]], py[t]],
-D[H[px[t], py[t], x[t], y[t]], x[t]],
-D[H[px[t], py[t], x[t], y[t]], y[t]]}

n = 140;
list = Range[n];
For[i = 1;, i < n + 1, i++,
ics = {Random[Real, {-0.55, 0.55}], Random[Real, {-0.45, 0.45}],
Random[Real, {-0.44, 0.9}], Random[Real, {-0.17, 0.9}]};
inisol = Solve[H[ini, ics[[4]], ics[[1]], ics[[2]]] == 0.866, ini];
inival = ini /. inisol;
ics = {ics[[1]], ics[[2]], inival[[2]], ics[[4]]};
TimeMin = 0;
TimeMax = 1000;
Section =
Reap[NDSolve[{x'[t] == a1, px'[t] == c1, y'[t] == b1, py'[t] == d1,
x[0] == ics[[1]], y[0] == ics[[2]], px[0] == ics[[3]],
py[0] == ics[[4]]}, {x[t], y[t], px[t], py[t]}, {t, TimeMin,
TimeMax}, MaxSteps -> \[Infinity],
Method -> {"EventLocator", "Event" -> x[t],
"EventAction" :> Sow[{y[t], py[t]}]}]][[2]];
list[[i]] =
ListPlot[Transpose[Table[Section, {j, Length[Section]}]][[1]][[1]],
PlotStyle -> {RGBColor[Random[], Random[], Random[]],
AbsolutePointSize[0.5]}]]

Show[list, PlotRange -> All, Frame -> True, Axes -> False,
ImageSize -> Large, LabelStyle -> Directive[White, Small],
Background -> GrayLevel[0.1]]


### Another Method

With this code, I get a different result

     With[{icv = {{0, 0.4082401254164024, 0, 0}, {0, 0.24008038662195985, -0.15316368171208566, 0.28838672839135326}, {0, 0.24352018790129148,         0.09848498536170303, 0.3135210500944904}, {0, 0.20788819727584018, 0.18687849877360435,
0.3047456328865395}, {0, 0.11124773284074371, 0.24234454959281115, 0.32410152798705066}, {0, 0.09612790055948518, 0.2562151235263743,        0.32091473672760257}, {0, 0.14504873906640348, 0.2433184249171519, 0.3098719090788399}, {0, 0.21466692008332666, 0.21341797098009158, 0.2855018063099425}, {0, 0.2242200636897092,  0.21057030870037158, 0.27976766448196655}, {0,        0.3683386516212832, 0.18778901046261137, 0.011696534024095057}, {0, 0.3451153373835715,  0.23762310035771178, 0.005962392196119196}, {0, 0.33504538978074816, 0.256132905175892, 0.0016617858251372995}, {0, 0.37097564355830726, -0.03575247849541091,        0.1665183633794433}, {0, 0.37507484064199426, -0.06280527015275122, 0.14788240243852177}, {0, 0.3774580162329546, -0.08416273725065143, 0.12924644149760023}, {0, 0.3779249463635693, -0.09982487978911159, 0.11491108692766056}, {0, 0.3818158505457779, -0.10979169776813169,
0.08910744870176919}, {0, 0.3850735819245613, -0.11406319118771174,
0.06617088138986574}, {0, 0.3571135001050445, -0.02151416709681081,
0.1966226079763166}, {0, 0.3537329768656702, -0.0015805311387706023,
0.20379028526128642}}},

psection =
Reap[Table[
NDSolve[{x'[t] == px[t],
px'[t] == -(2 g^2 v^2 x[t] + g^2 x[t] y[t]^2), y'[t] == py[t],
py'[t] == -(2 g^2 v^2 y[t] + g^2 x[t]^2 y[t]),
x[0] == Part[icv[[i]], 1], px[0] == Part[icv[[i]], 2],
y[0] == Part[icv[[i]], 3], py[0] == Part[icv[[i]], 4]}, {x,
px, y, py}, {t, 0, 1000}, MaxSteps -> \[Infinity],
Method -> {"EventLocator", "Event" -> x[t],
"EventAction" :> Sow[{y[t], py[t]}]}], {i, 1,
Length[icv]}]][[2]];]

ListPlot[psection, PlotRange -> All, PlotTheme -> "Scientific",
LabelStyle -> Directive[Black, 18],
FrameLabel -> {{"\!$$\*SubscriptBox[\(p$$, $$y$$]\)", None}, {"y",
"YMH system"}}, RotateLabel -> False,
Epilog -> {Text["E=10", {0.3, 0.4}]}, ImageSize -> 400]


I think there is a problem in defining the energy surface.

• I cannot find the correct plots you show in any of the two papers you link ... Commented Sep 14, 2023 at 14:46
• @Domen .. Please see this link sci-hub.hkvisa.net/10.1142/S0217732397001503 or this doi.org/10.1142/S0217979299002447 Commented Sep 14, 2023 at 15:00
• I've edited the question. Commented Sep 14, 2023 at 15:04
• (1) Look at the range of $q_2$ on the original plots, it spans $-3\lesssim q_2 \lesssim 3$. Yet, you sample $q_2$ only from {-0.45, 0.45}. Increase this range (also for $p_2$). (2) You look at $E=0.866$ (why?), when the original plots are for $E=10$. After fixing (1) and (2), you should get the correct image. Commented Sep 14, 2023 at 15:25
• Could you share your code, please? I'm getting imaginary numbers when I do the necessary corrections as you suggest. This was the initial reason for taking them as such Commented Sep 14, 2023 at 16:28

I'm currently co-developing an interactive tool for extracting Poincaré sections from Hamiltonian systems. I've got a preliminary version up and running, so feel free to use it. It's available on the Wolfram Function Repository as ClickPoincarePlot2D. For your particular problem, the initial Poincaré section with $$v=1$$, $$g=1$$, and energy $$E=10$$ would be as follows:

H := Function[{S}, 1/2 (S[[2]]^2 + S[[4]]^2) + S[[1]]^2 + S[[3]]^2 + 1/2*S[[1]]^2*S[[3]]^2]

hameqns = {x'[t] == px[t], px'[t] == -(2 x[t] + x[t]*y[t]^2),
y'[t] == py[t], py'[t] == -(2 y[t] + x[t]^2*y[t])};

style = {PlotStyle -> {{AbsolutePointSize[1], GrayLevel[0], Opacity[0.4]}},
AspectRatio -> 1, PlotHighlighting -> None};
(*Omit the PlotHighlighting option if you're using a version prior to 13.3.*)

ResourceFunction["ClickPoincarePlot2D"][hameqns, H, 10, t, 1000, x[t], {y[t], py[t]}, style]


The output with the Poincaré section looks as follows:

I'll be continuing the development of this tool to address its limitations and will be adding features such as zooming into specific chaotic regions, allowing users to click on stability islands and display additional orbits. The advantage of this tool is its reduced resource consumption compared to the other code I provided in response to a different question about Poincaré sections.

You can supplement this tool with the Wolfram Community post titled Plot Poincare Section of a given Hamiltonian:

PoincareSections[eqns_?VectorQ, depvars_?VectorQ, cross_,
retrievesec_?VectorQ, icv : {__?NumericQ}, tmax_Integer?Positive] /;
Length[retrievesec] == 2 && Length[icv] == 4 :=
Module[{vars, time},
time = Times @@ Variables@Level[depvars, {-1}];
Reap[NDSolve[
Join @@ {eqns,
Sow[retrievesec, "LocationMethod" -> "EvenLocator"]]}},
vars, {time, 0, tmax}, MaxSteps -> ∞,
MaxStepSize -> 0.1]][[-1, 1]]]

PoincareSections[eqns_?VectorQ, depvars_?VectorQ, cross_,
retrievesec_?VectorQ, icv_?MatrixQ, tmax_Integer?Positive] :=
Map[PoincareSections[eqns, depvars, cross, retrievesec, #, tmax] &, icv]


Next, I provide the initial conditions obtained with the interactive routine:

    icvals = {{0, 4.462528320177665, 0.17303856765541714,
0.16110958947199644}, {0, 4.467736758629249,
0.06879141200183614, 0.1728112810629193}, {0,
4.4629142781672035, 0.08539567409649612,
0.2604060356875785}, {0, 4.457078219133263,
0.09093042812804947, 0.34339054006883457}, {0,
4.44273956486575, 0.11306944425426274,
0.48630829761433114}, {0, 4.432274545806214,
0.1094320279123828, 0.5753187057034409}, {0, 4.422434117182368,
0.10295229672798317, 0.6487512068414824}, {0,
4.4126874534322456, 0.11523870068793313,
0.7082580899990921}, {0, 4.401848022481288,
0.12845252714010807, 0.768592143805943}, {0, 4.388180737932256,
0.14463303584951864, 0.8378737262393601}, {0,
4.372830596551317, 0.16618755584494688,
0.9072573870899112}, {0, 4.354667994170191,
0.18369170535991586, 0.9845716709888115}, {0,
4.335255458734564, 0.20905045373732956, 1.057428921061173}, {0,
4.381078248441183, 0.4893283115828355,
0.5720742853882897}, {0, 4.311513492521181, 0.2385748806574618,
1.1388658641149816}, {0, 4.280339885860469,
0.09636574678509191, 1.2884555666398008}, {0,
4.255797192593596, 0.09053334214390507,
1.3681365733847572}, {0, 4.283226359375878, -0.46353021338614,
1.1064589630456512}, {0, 4.3331024584312, -0.6526232829598253,
0.6102366637747763}, {0,
4.042356859930489, -0.09604607887886041,
1.9081145978260585}, {0, 2.122007161924683,
0.08359063563407387, 3.934858423883753}, {0,
3.6138627754729127, -1.2673396781350363,
1.9307242476595103}, {0,
3.5047138061488123, -1.3987836406681824,
1.9503306873933028}, {0,
3.263100249458284, -1.3390363849712978,
2.4012788012705393}, {0,
3.3213488705931584, -1.2075924224381516,
2.460098120471918}, {0, 3.460895201964421, -0.9925023019293671,
2.460098120471918}, {0,
2.874602524936952, -1.2075924224381516, 2.969865553550533}, {0,
3.038063255679484, -1.016401204208121, 2.95025911381674}, {0,
3.5631421293932606, 0.47728018821399454,
2.6169496383422612}, {0,
3.1217459031207135, -0.6937660234449439,
3.0482913124857043}, {0,
3.1145655514289583, -0.27553523356675136,
3.1855363906222545}, {0,
1.449807234660935, -1.9843067464976518, 3.165929950888462}, {0,
2.0900579174201654, -2.79568021014104, \
-0.0014923407056293622}, {0,
1.832519509623928, -2.8844681909476537,
0.03948110522791894}, {0, 1.503596035819526, 2.976819841668624,
0.12761748141920554}, {0,
3.5415459321081486, -0.8960177933209614,
2.4190404376266006}, {0,
3.130688321590937, -0.07791485738605086,
3.191653048664773}, {0, 3.5390581462779105, 0.8440741656517057,
2.4597042592601888}, {0,
2.248738996441609, -2.5581951869664943, 1.361781075153312}, {0,
2.4662518446386414,
2.5711803637365174, -0.8340652930604416}, {0,
2.4242287391373676,
2.558194602848661, -1.0170524904115879}, {0, 2.110487349745546,
2.5711803637365174, 1.524436361687664}, {0, 2.351295468570382,
2.5711803637365174, 1.1177981453517838}, {0,
2.6184146448394854,
2.519237320185094, -0.6714100065260895}, {0,
2.8320631190315435, 2.4413227548579597,
0.24352598022964114}, {0, 4.177374473484878,
0.2856864474739096, 1.544768272504458}, {0,
4.172899517751086, -0.064929096498195, 1.60576400495484}, {0,
3.8439427761413656, -0.9869181195359515, -1.8099970122665545}, \
{0, 3.927220540238002, -1.1167757284145088, 1.443108718420488}, {0,
2.805662794441176, 1.7141201451380388, 2.500368080893777}, {0,
3.513796576017365, -1.4673912723866134,
1.8294150239395743}, {0, 3.160908672942136, -1.207676054629499,
2.6630233674281274}, {0,
4.290890526341124, -0.8570605106573943,
0.34518553431360965}, {0,
4.241737814626106, -0.14284366182532926,
1.4024448967868983}, {0, 1.2944984720466737,
3.0256819948114675, 0.12153451532887552}, {0,
0.9972003079434875, 3.0821237104495274,
0.0813536281126289}, {0,
1.924047846028568, -2.8439473819482592,
0.3492369734322397}, {0,
2.240082902959127, -2.7366286067454566,
0.06128886811690352}, {0, 2.1940388127188575,
2.750043775497823, 0.24639836439104817}};


Now we use the PoincareSections routine with the above initial conditions:

hameqns = {x'[t] == px[t], px'[t] == -(2 x[t] + x[t]*y[t]^2),
y'[t] == py[t], py'[t] == -(2 y[t] + x[t]^2*y[t])};
states = {x[t], px[t], y[t], py[t]};
cross = x[t];
retrievesec = {y[t], py[t]};

ListPlot[PoincareSections[hameqns, states, cross, retrievesec, icvals, 3000],
PlotStyle -> {AbsolutePointSize[0.05]}, Axes -> False,
Frame -> True, FrameStyle -> Directive[White, Medium],
Background -> GrayLevel[0], ImageSize -> Full]


The Poincaré section of Yang-Mills-Higgs System for $$v=1$$, $$g=1$$, and $$E=10$$:

Another approach using ClickPane:

YMHSection[{x0_, px0_, y0_, py0_}, {v_, g_}, e_, tmax_] :=
Block[{dx0, x, px, y, py},
dx0 = Sqrt[2*e - py0^2 - 2*g^2*v^2*x0^2 - 2*g^2*v^2*y0^2 - g^2*x0^2*y0^2];
If[# == {}, {}, First[#]] &@Last[Reap[NDSolve[{
x'[t] == px[t],
px'[t] == -(2 g^2*v^2* x[t] + g^2 x[t]* y[t]^2),
y'[t] == py[t],
py'[t] == -(2 g^2*v^2*y[t] + g^2 x[t]^2*y[t]),
x[0] == x0,
px[0] == dx0,
y[0] == y0,
py[0] == py0,
WhenEvent[x[t] == 0,
Sow[{y[t], py[t]}, "LocationMethod" -> "EvenLocator"]]
},
{x, px, y, py}, {t, 0, tmax}, MaxSteps -> ∞,
MaxStepSize -> 0.1
]]]]]


Poincaré sections considering $$v=1$$, $$g=1$$ and $$E=10$$:

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,
YMHSection[{0, 0, #[[1]], #[[2]]}, {1, 1}, 10, 2000]]) &],
Button[Style["Copy orbit list", 12],
CopyToClipboard[Iconize[f, "Orbit list"]],
ImageSize -> Automatic, Alignment -> Center,
Appearance -> {"Palette", "Pressed"}]}, "  "]},
Alignment -> Center]]


I hope this is helpful and that you can plot some excellent Poincaré sections.

• It's wonderful. Thank you for sharing. However, I'm using Mathematica 12.0 and it has some compatibility issues, I guess. Error : "PlotHighlighting::shdw: Symbol PlotHighlighting appears in multiple contexts {System,Global}; definitions in context System may shadow or be shadowed by other definitions." Commented Sep 14, 2023 at 16:41
• Don't worry, I'll set up a more hands-on interactive version for your issue using ClickPane.ClickPane is also available in the version you're using. Commented Sep 14, 2023 at 16:43
• I forgot to mention that you can omit the PlotHighlighting option, as it's not available in your version. This way, you can use the ClickPoincarePlot2D resource function without issues and the other version using ClickPane. See the update, please! :-) Commented Sep 15, 2023 at 1:14
• @E.Chan-López It is a nice tool (+1). How you compute "displaying initial values" at fixed Hamiltonian? In your example there are only zeroes in the first column. Commented Sep 15, 2023 at 3:51
• @E.Chan-López Please see update to my post how to compute initial data with using NMinimize. Commented Sep 15, 2023 at 13:16

If we take icv from your second example and define Hamiltonian as

H0[{x_, px_, y_, py_}] :=
1/2 (px^2 + py^2) + g^2*v^2*(x^2 + y^2) + 1/2*g^2*x^2*y^2;


then we have for energy levels

Map[H0, icv]

(*{0.08333, 0.0938619, 0.0884981, 0.102967, 0.11744, 0.12176, 0.117734, \
0.109344, 0.108612, 0.10317, 0.116035, 0.121733, 0.0839539, \
0.0852197, 0.086673, 0.0879809, 0.088916, 0.0893405, 0.0835581, \
0.0833312}*)


Therefore in the first example you used $$H=0.866$$ while in the second example $$0.08333\le H \le 0.12176$$. Also initial data in the first example are not same as in the second. Please, note, that in the second example we can use icv and -icv as well as follows

g = 1; v = 1.0(*1.1*)(*1.2*); H =
1/2 (px[t]^2 + py[t]^2) + g^2*v^2*(x[t]^2 + y[t]^2) +
1/2*g^2*x[t]^2*y[t]^2; {a1, b1, c1, d1} = {D[H, px[t]],
D[H, py[t]], -D[H, x[t]], -D[H, 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 + 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]];

icv = {{0, 0.4082401254164024, 0, 0}, {0,
0.24008038662195985, -0.15316368171208566,
0.28838672839135326}, {0, 0.24352018790129148,
0.09848498536170303, 0.3135210500944904}, {0,
0.20788819727584018, 0.18687849877360435,
0.3047456328865395}, {0, 0.11124773284074371,
0.24234454959281115, 0.32410152798705066}, {0,
0.09612790055948518, 0.2562151235263743,
0.32091473672760257}, {0, 0.14504873906640348,
0.2433184249171519, 0.3098719090788399}, {0,
0.21466692008332666, 0.21341797098009158,
0.2855018063099425}, {0, 0.2242200636897092,
0.21057030870037158, 0.27976766448196655}, {0,
0.3683386516212832, 0.18778901046261137,
0.011696534024095057}, {0, 0.3451153373835715,
0.23762310035771178, 0.005962392196119196}, {0,
0.33504538978074816, 0.256132905175892,
0.0016617858251372995}, {0,
0.37097564355830726, -0.03575247849541091,
0.1665183633794433}, {0,
0.37507484064199426, -0.06280527015275122,
0.14788240243852177}, {0,
0.3774580162329546, -0.08416273725065143,
0.12924644149760023}, {0,
0.3779249463635693, -0.09982487978911159,
0.11491108692766056}, {0,
0.3818158505457779, -0.10979169776813169,
0.08910744870176919}, {0,
0.3850735819245613, -0.11406319118771174,
0.06617088138986574}, {0,
0.3571135001050445, -0.02151416709681081,
0.1966226079763166}, {0,
0.3537329768656702, -0.0015805311387706023,
0.20379028526128642}}; iniv = Join[icv, -icv];

abcdata = Map[psect, iniv];


Visualization

ListPlot[abcdata, ImageSize -> Medium, Background -> GrayLevel[0.1],
Frame -> True]


Update 1.To compute initial data at fixed Hamiltonian H=E0 we can use NMinimize as follows

g = v = 1;
h[x_, px_, y_, py_] :=
1/2 (px^2 + py^2) + g^2*v^2*(x^2 + y^2) + 1/2*g^2*x^2*y^2;

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[(h[x, px[[i]], y, py] - E0)^2, {x, py, y}][[2]];, {i,
n}]; Transpose[{X, px, Y, P}]];


In a case of H=10, n=20 we have, for example,

iniv = inivalues[10, 20, -.5, .5]


Mapping Hamiltonian in the form H0 we can check that

Map[H0, iniv]

(*Out[]= {10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10.}*)


Finally we compute section data and plot

abcdata = Map[psect, iniv];

ListPlot[abcdata, ImageSize -> Large, Background -> GrayLevel[0.05],
Frame -> True]


• (+1) This is excellent, @AlexTrounev! Utilizing NMinimize is the most intelligent approach for computing the initial conditions constrained by the energy level restriction. This mitigates numerous issues within a non-interactive routine and is highly likely to, when suitably adapted to an interactive routine, address some of the challenges I observe with NSolve in more complex systems. Commented Sep 15, 2023 at 14:52
• @E.Chan-López Thank you. Please note, that NMinimize is very effective in a case of convex functions like Hamiltonian discussed above. Commented Sep 15, 2023 at 15:29
• What you did gives me ideas to improve the routine. Thanks, @AlexTrounev! I'll happily cite this post and what you did in the next version of the routine. I hope collaboration is possible. Commented Sep 15, 2023 at 15:58
• Thank you, I'm glad this is useful for you. Commented Sep 16, 2023 at 5:20
• I managed to create a compatible structure for the approach you shared with me where you use NMinimize. I'm thinking about integrating it into the interactive routine. Would you be interested in collaborating on an article about this? We could document the feature and perhaps delve into a couple more complex examples. Additionally, I have a zoom function to click and draw orbits in the magnified region. With these additions, our routine could be super solid. Commented Oct 4, 2023 at 15:19