Plot of massless particles on the unstable null circular geodesic

Question

My friend and I were trying to reproduce Figure 4(left), page 7, from the paper "Probing phase structure of black holes with Lyapunov exponents", Guo 2022. We succeeded in reproducing the result after writing a lengthy and complicated code. Here I dissect the code in steps:

Step 1:
Find the parametric plot for Free Energy v/s Temperature and the Transition Temperature.

l = 1;
T[rp_, Q_] := 1/(4 \[Pi] rp) (1 - Q^2/rp^2 + (3 rp^2)/l^2)
M[rp_, Q_] := rp/2 (1 + Q^2/rp^2 + rp^2/l^2)
S[rp_] := \[Pi] rp^2
\[Phi][rp_, Q_] := Q/rp
F[rp_, Q_] := Simplify[M[rp, Q] - (T[rp, Q] S[rp])]
Q = Qt l;
rp = rpt l;
Tt[rpt_, Qt_] = T[rp, Q] l;
Ft[rpt_, Qt_] = F[rp, Q]/l;
Mt[rpt_, Qt_] = M[rp, Q]/l;
r = rt l;
pplt[Qt_] := 
 ParametricPlot[Evaluate@{Tt[rpt, Qt], Ft[rpt, Qt]}, {rpt, 0.01, 3}, 
 PlotRange -> {{0.2, 0.6}, {-0.05, 0.15}}, AspectRatio -> 1]
ip[Qt_] := Graphics`Mesh`FindIntersections[pplt[Qt]][[1]]
Tp = Table[ip[Qt][[1]], {Qt, 0.01, 0.166, 0.001}];

Step 2:
Find the critical values and generate a Table for Qtt(Charge).

rptc = SolveValues[D[Tt[rpt, Qt], {rpt, 2}] == 0, rpt][[2]];
Qtc = SolveValues[D[Tt[rpt, Qt], rpt] == 0, Qt][[2]];
rpt1 = rptc /. Qt -> Qtc;
rptc = Last @@ SolveValues[rpt1 == rpt, rpt];
Qtc = Qtc /. rpt -> rptc;
Ttc = Tt[rpt, Qt] /. rpt -> rptc /. Qt -> Qtc;
\[Phi]c = \[Phi][rp, Q] /. Qt -> Qtc /. rpt -> rptc;
f[r_, rpt_, Qt_] := 1 - (2 Mt[rpt, Qt])/r + Q^2/r^2 + r^2/l^2
Veff[r_] := f[r, rpt, Qt] (L^2/r^2 + \[Delta]1)
r0 = l/2 (3/2 rpt (Qt^2/rpt^2 + rpt^2 + 1) + 
 Sqrt[9/4 rpt^2 (Qt^2/rpt^2 + rpt^2 + 1)^2 - 8 Qt^2]);
l = 1; \[Delta]1 = 0;
\[Lambda][rpt_, Qt_] = 
 Sqrt[-((r0^2 f[r0, rpt, Qt])/L^2) Veff''[r0]] // Simplify;
Qtt = Table[Qt, {Qt, 0.01, 0.166, 0.001}];

Step 3:
Find the parametric plot for Lyapunov exponent v/s Temperature and find the intersection points.

pplt1[Qt_] := 
ParametricPlot[{Tt[rpt, Qt], \[Lambda][rpt, Qt]}, {rpt, 0.01, 3}, 
PlotRange -> {{0.0, 0.6}, {0.75, 26.00}}, AspectRatio -> 1, 
AxesLabel -> {"T", "\[Lambda]"}, PlotPoints -> 500]
pts = Table[
 Graphics`Mesh`FindIntersections[
 Show[pplt1[Qtt[[t]]], 
 Graphics[{Red, Line[{{Tp[[t]], 1}, {Tp[[t]], 26}}]}]]], {t, 1, 
 Length[Tp]}];

which gives the following result after using animate: moving line
Step 4:
Find the minimum and maximum intersection points and find their difference to get the raw plot:

minmax = ReplacePart[pts, {{_, _, 1} -> Nothing, {_, 2 | 2} -> Nothing}];
\[CapitalDelta]\[Lambda] = Flatten[Differences[#] & /@ minmax];
data = Transpose@{Tp, \[CapitalDelta]\[Lambda]};
ListLinePlot[data, PlotRange -> All]

Step 5
Find the rescaled(normalized) plot using the critical values.

\[Lambda]c = \[Lambda][rptc, Qtc] // N;
diff\[Lambda] = \[CapitalDelta]\[Lambda]/\[Lambda]c;
t1 = Tp/Ttc;
data2 = Transpose@{t1, diff\[Lambda]};
ListLinePlot[data2, PlotRange -> All]

I wish to use functional programming techniques to reduce the code length to get the desired result. Is there a way to shorten the code? Any help in this regard would be truly beneficial!

• 200_success: What do you mean by "I wish to automate the process"?
• codebpr: @200_success by that I meant that I would like to remove some intermediate steps. I have also edited the question to directly address my specific concern.
• Sᴀᴍ Onᴇᴌᴀ: I changed the title so that it describes what the code does per site goals: "State what your code does in your title, not your main concerns about it.". Feel free to edit and give it a different title if there is something more appropriate.