Definitions

SetDirectory[NotebookDirectory[]];

SetAttributes[untilSuccess,HoldAll];
untilSuccess[e_]:=Quiet[Check[e,untilSuccess[e]]]

μm[f_]:=(Sqrt[4f''[1]])

Eopt[f_,μ_]:=((μ(f[1] f'[1]-((f'[1])^2)-f[1] f''[1]))/(2 f'[1] f''[1]))

Eth[f_]:=Eopt[f,μm[f]]

Eμopt[f_,μ_]:=-((μ f'[1])/(f'[1]+f''[1]))

𝒟[f_][μ_]:=With[{μm=(Sqrt[4f''[1]])},(1/2) +(((μ )^2)/((μm)^2))+Log[(1/2)μm]-UnitStep[(μ^2)-((μm)^2)](Log[RealAbs[(μ/(μm))]-(Sqrt[((μ^2)/((μm)^2))-1])]+ RealAbs[(μ/(μm))](Sqrt[((μ^2)/((μm)^2))-1]))]
d𝒟[f_][μ_]:=With[{μm=(Sqrt[4f''[1]])},(2/(μm))((μ/(μm))-UnitStep[(μ^2)-((μm)^2)]RealSign[μ] (Sqrt[((μ^2)/((μm)^2))-1]))]

Σa[f_][E_,μ_]:=𝒟[f][μ]-(1/2)-(1/2)Log[f'[1]]-(1/(2f[1]))((E^2)+((((f[1] μ+f'[1] E))^2)/(f[1] (f'[1]+f''[1])-((f'[1])^2))))

ΣaSaddles[f_][EE_,μ_]=Simplify[Σa[f][EE,μ],Assumptions->{(μ^2)<((μm[f])^2),μ>0}];
ΣaMinima[f_][EE_,μ_]=Simplify[Σa[f][EE,μ],Assumptions->{(μ^2)>((μm[f])^2),μ>0,f''[1]>0}];

𝒮tmp1[f_][E_,μ_][q1_,x_][β_,rd_,r1_,dd_,d1_]:=With[{Δx=1-x},
β E-rd μ+(1/2)((β^2)(f[1]-Δx f[q1])+(2 β rd-dd)f'[1]-Δx(2 β r1-d1)f'[q1]+(rd^2)f''[1]-Δx (r1^2)f''[q1]+(1/x)Log[(((rd-Δx r1))^2)+dd(1-Δx q1)-Δx d1(1-Δx q1)]-((Δx)/x)Log[(((rd-r1))^2)+(dd-d1)(1-q1)])
]

𝒮tmp2[f_][EE_,μ_][q1_,x_][β_,rd_,r1_]=Simplify[𝒮tmp1[f][EE,μ][q1,x][β,rd,r1,dd,d1]/.First[Solve[0==D[𝒮tmp1[f][EE,μ][q1,x][β,rd,r1,dd,d1],{{dd,d1}}],{dd,d1}]],Assumptions->{((1-q1 (1-x))/(f'[1]-(1-x) f'[q1]))>0,((1-q1)/(f'[1]-f'[q1]))>0,0<x<1,0<q1<1}];

𝒮[f_][EE_,μ_][q1_,x_]=𝒟[f][μ]+Simplify[𝒮tmp2[f][EE,μ][q1,x][β,rd,r1]/.First[Solve[0==D[𝒮tmp2[f][EE,μ][q1,x][β,rd,r1],{{β,rd,r1}}],{β,rd,r1}]]];

u[f_]:=f[1](f'[1]+f''[1])-((f'[1])^2)
v[f_]:=f'[1](f''[1]+f'''[1])-((f''[1])^2)
w[f_]:=2f''[1](f''[1]-f'[1])+f'[1]f'''[1]
y[f_]:=f'[1](f'[1]-f[1])+f''[1]f[1]
z[f_]:=f[1](f''[1]-f'[1])+((f'[1])^2)

μ0[f_][E_]:=((-1)/(z[f]))(2E f'[1]f''[1]+(Sqrt[2f''[1] u[f](Log[((f''[1])/(f'[1]))]z[f]-(E^2)(f''[1]-f'[1]))]))

Emax[f_]=EE/.Solve[μ0[f][EE]==μm[f],EE][[1]]//Simplify

q1Exp={f[q1]->f[1]+f'[1](q1-1)+(1/2)f''[1](((q1-1))^2),Derivative[n_][f][q1]:>Derivative[n][f][1]+Derivative[n+1][f][1](q1-1)+(1/2)Derivative[n+2][f][1](((q1-1))^2)};

d𝒮dq1dx=Simplify[Limit[Simplify[D[D[𝒮[f][ϵ,μ][q1,x],q1],x]/.x->1]/.q1Exp,q1->1,Direction->"FromBelow"]];

d𝒮Coeff=FullSimplify[CoefficientList[d𝒮dq1dx/.μ->((-1)/(z[f]))(2ϵ f'[1]f''[1]+2y),{ϵ,y}]]

ΣaCoeff=FullSimplify[CoefficientList[(2 f[1] f'[1]-2 ((f'[1])^2)-2 f[1] f''[1])ΣaSaddles[f][ϵ,((-1)/(z[f]))(2ϵ f'[1]f''[1]+2y)],{ϵ,y}]]

letterRules={d->4d𝒮Coeff[[1,1]],c->2d𝒮Coeff[[2,2]],b->4d𝒮Coeff[[3,1]],a->4d𝒮Coeff[[1,3]],g->-ΣaCoeff[[1,1]],e->ΣaCoeff[[3,1]],h->ΣaCoeff[[1,3]]/2};

solLetters=Solve[{0==a (y^2)+b (ϵ^2)+2c ϵ y+d,0==2h (y^2)+e (ϵ^2)-g},{y,ϵ}]//Simplify;

solLetters//FullSimplify

letterRules

E1RSBp[f_]:=Module[{
a=(-6 ((y[f])^2) (f''[1]-f'[1]) f''[1]+w[f] (3 ((y[f])^2)-4f[1]f'[1]f''[1](f'[1]-f[1])))/( f'[1](((f''[1]u[f]z[f]))^2)),
b=((f'[1] w[f])/((z[f])^2)),
c=((w[f])/(f''[1] ((z[f])^2))),
d=(1/(f'[1]f''[1]))w[f],
e=f''[1]-f'[1],
g=Log[f''[1]/f'[1]] z[f],
h=(1/(f''[1] u[f]))
},
{-(c g+(Sqrt[(c^2) (g^2)+(2d h-a g) (b g-d e)]))(((2 (c^2) e g+(2 b h-a e) (b g-d e)+2 c e(Sqrt[(c^2) (g^2)+(2d h-a g) (b g-d e)])))^(-1/2)) ,
-Sign[Log[((f''[1])/(f'[1]))] ((f'[1])/(z[f]))-((1/(f'[1]))-(1/(f''[1])))](c  g-(Sqrt[(c^2) (g^2)+(2d h-a g) (b g-d e)]))(((2 (c^2) e g+(2 b h-a e) (b g-d e)-2 c e(Sqrt[(c^2) (g^2)+(2d h-a g) (b g-d e)])))^(-1/2))}
]

G[f_]:=f'[1]Log[((f''[1])/(f'[1]))](3y[f](f''[1]-f'[1])f'''[1]-2(f'[1]-2f[1])f''[1]w[f])-2(f''[1]-f'[1])u[f]w[f]-2((Log[((f''[1])/(f'[1]))])^2)((f'[1])^2)f''[1]v[f]

h11Auff[f_,x_]:=((x f'[1]-f'[x])/(f'[1]-f'[x]))+(((((f'[1]x-f'[x]))^2)f[x])/(f'[1]f'[x](f'[1]-f'[x])x(1-x)))-Log[((f'[1]-f'[x])/(f'[1](1-x)))]
h12Auff[f_,x_]:=(1-x)(f'[1]-f'[x])(((f'[1]f'[x]x)/(((f'[1]x-f'[x]))^2))Log[((f'[1]x)/(f'[x]))]-((f'[x])/(f'[1]x-f'[x])))+f'[x](1-x)-1+f[x]
Ψ[λ_,p_,s_]:=With[{f=Function[q,(λ (q^p)+(1-λ)(q^s))]},-((f'[1]-1)/(f''[1]))-Log[((f''[1])/(f'[1]))]+1-(2/(f'[1]))+((f''[1])/((f'[1])^2))]

cleanContourPlot[cp_Graphics]:=Module[{points,groups,regions,lines},groups=Cases[cp,{style__,g_GraphicsGroup}:>{{style},g},Infinity];
points=First@Cases[cp,GraphicsComplex[pts_,___]:>pts,Infinity];
regions=Table[Module[{group,style,polys,edges,cover,graph},{style,group}=g;
polys=Join@@Cases[group,Polygon[pt_,___]:>pt,Infinity];
edges=Join@@(Partition[#,2,1,1]&/@polys);
cover=Cases[Tally[Sort/@edges],{e_,1}:>e];
graph=Graph[UndirectedEdge@@@cover];
{Sequence@@style,FilledCurve[List/@Line/@First/@Map[First,FindEulerianCycle/@(Subgraph[graph,#]&)/@ConnectedComponents[graph],{3}]]}],{g,groups}];
lines=Cases[cp,_Tooltip,Infinity];
Graphics[GraphicsComplex[points,{regions,lines}],Sequence@@Options[cp]]]

cleanContourPlot[Legended[cp_Graphics,rest___]]:=Legended[cleanContourPlot[cp],rest]
cleanRegionPlot[cp_Graphics]:=Module[{points,groups,regions,lines},groups=Cases[cp,{style__,g_GraphicsGroup}:>{{style},g},Infinity];
points=First@Cases[cp,GraphicsComplex[pts_,___]:>pts,Infinity];
regions=Table[Module[{group,style,polys,edges,cover,graph},{style,group}=g;
polys=Join@@Cases[group,Polygon[pt_,___]:>pt,Infinity];
edges=Join@@(Partition[#,2,1,1]&/@polys);
cover=Cases[Tally[Sort/@edges],{e_,1}:>e];
graph=Graph[UndirectedEdge@@@cover];
{Sequence@@style,FilledCurve[List/@Line/@First/@Map[First,FindEulerianCycle/@(Subgraph[graph,#]&)/@ConnectedComponents[graph],{3}]]}],{g,groups}];
lines=Cases[cp,{__,_Line},Infinity];(*only change*)Graphics[GraphicsComplex[points,{regions,lines}],Sequence@@Options[cp]]];

eqZero[f_,k_]:=-(1/2)Rd (OverHat[β])f'[1]-(1/2)(((OverHat[β])^2)(f[1]+Sum[(x[i]-x[i+1])f[q[i]],{i,0,k}])+(1/(x[1]))Log[1+(OverHat[β])(Rd^(-1))(1+Sum[(x[i]-x[i+1]) q[i],{i,0,k}])]+Sum[((1/(x[j+1]))-(1/(x[j])))Log[1+(OverHat[β])(Rd^(-1))(1+Sum[(x[i]-x[i+1]) q[i],{i,j+1,k}]-x[j+1]q[j])],{j,1,k}])/.{x[0]->0,q[0]->0,x[k+1]->1}

nMaximize[f_,k_,n_]:=
First@Reverse@SortBy[
ParallelTable[untilSuccess@FindMaximum[Prepend[Flatten@Join[{Rd>0,(OverHat[β])>0,0<q[1],0<x[1],q[k]<1,x[k]<1},Table[{0<q[i+1]-q[i],0<x[i+1]-x[i]},{i,1,k-1}]],(1/(OverHat[β]))eqZero[f,k]],Thread[{Join[x/@Range[k],q/@Range[k],{(OverHat[β]),Rd}],RandomReal[{0,1},2k+2,WorkingPrecision->50]}],WorkingPrecision->50],{n}],
First]

fps[p_,s_][λ_]:=Function[q,(1/2)(λ (q^p)+(1-λ)(q^s))]

Ealg[f_]:=-Integrate[((f''[q])^(1/2)),{q,0,1}]

Ealg3s[s_]=-Integrate[((fps[3,s][λ]''[q])^(1/2)),{q,0,1},Assumptions->{s>3,λ>0,λ<1}]

FindRoot[{Emax[fps[3,s][λ]]-E1RSBp[fps[3,s][λ]][[1]],D[Emax[fps[3,s][λ]]-E1RSBp[fps[3,s][λ]][[1]],λ]/.Sign'[_]:>0},{{λ,0.5,(10^(-2)),1-(10^(-2))},{s,10}},WorkingPrecision->50]

Figures

Phases figure

λAuffStart=FindRoot[Evaluate[{Ψ[λ,3,s],D[Ψ[λ,3,s],λ]}],{{s,12.2},{λ,0.86}},WorkingPrecision->50]

λ12=Prepend[Rest[FoldList[Function[{sol,ss},Prepend[FindRoot[{h11Auff[Function[q,(λ (q^3)+(1-λ)(q^ss))],x],h12Auff[Function[q,(λ (q^3)+(1-λ)(q^ss))],x]},Rest[sol]/.Rule->List,MaxIterations->1000,WorkingPrecision->20],s->ss]],{s->(s/.λAuffStart),λ->(0.8901),x->0.9664},s+Rest@Range[0,8,1/100]/.λAuffStart]],λAuffStart];

λ21=Prepend[Rest[FoldList[Function[{sol,ss},Prepend[FindRoot[Ψ[λ,3,ss],Rest[sol]/.Rule->List,MaxIterations->1000,WorkingPrecision->10],s->ss]],{s->(s/.λAuffStart),λ->(λ+0.01/.λAuffStart)},s+Rest@Range[0,8,1/100]/.λAuffStart]],λAuffStart];

rsbComplexityStart=FindRoot[D[G[fps[3,s][λ]],λ]/.λ->0,{s,4.5},WorkingPrecision->20]

rsbComplexityBoundary=Prepend[Rest[FoldList[Function[{sol,ss},Prepend[FindRoot[G[Function[q,(1/2)(λ (q^3)+(1-λ)(q^ss))]],Rest[sol]/.Rule->List,WorkingPrecision->15],s->ss]],Append[rsbComplexityStart,λ->(10^(-1))],s+Rest[Range[0,16,1/100]]/.rsbComplexityStart]],Append[rsbComplexityStart,λ->0]];

plot3s=Show[
ListPlot[{{s,λ}/.rsbComplexityBoundary,{s,λ}/.λ12,{s,λ}/.λ21,{{8,0},{8,1.1},{20,1.1}}},Joined->True,Filling->{1->0,2->{3},4->0},PlotStyle->{ColorData[97][1],ColorData[97][3],ColorData[97][3],ColorData[97][2]},AspectRatio->1/GoldenRatio,PlotRange->{{3,17},{0,1}},FrameLabel->{s,λ},RotateLabel->False,LabelStyle->Directive[Black,FontFamily->"Times",FontSize->10,Filling->{{1}->0}],FrameStyle->Black,Frame->True,FillingStyle->Opacity[0.1]],
Epilog->{Text[Style["RSB",Directive[ColorData[97][1],FontFamily->"Times",FontSize->10]],{6.32,0.4+0.035}],
Text[Style["complexity",Directive[ColorData[97][1],FontFamily->"Times",FontSize->10]],{6.32,0.4-0.035}],
Text[Style["Non\[Hyphen]convex χ",Directive[ColorData[97][2],FontFamily->"Times",FontSize->10]],{12.5,0.4}],
Text[Style["RSB T = 0",Directive[ColorData[97][3],FontFamily->"Times",FontSize->10]],{15,0.85}]},ImageSize->270,RotateLabel->False
]

Export["figs/phases_34.pdf",plot3s];

Complexity figure

f35=Function[q,(1/2)((1/2) (q^3)+(1-(1/2))(q^5))];

Eini=E1RSBp[f35][[2]]-(10^(-3));

zeroComplexityEqs=Append[D[𝒮[f35][ϵ,μ][q1,x],{{q1,x}}],𝒮[f35][ϵ,μ][q1,x]];

μ0[f35][Eini]//N

iniComplexityZero=Select[Table[FindRoot[zeroComplexityEqs/.ϵ->Eini,{{μ,μ0[f35][Eini]},{q1,1-0.01RandomReal[],0,1},{x,1-0.01RandomReal[],0,1}},WorkingPrecision->15],{100}],Norm[zeroComplexityEqs/.ϵ->Eini]<(10^(-4))&&x<1/.#&];

Length[iniComplexityZero]

iniComplexityZero2=FindRoot[zeroComplexityEqs/.ϵ->Eini,(#1->Re[#2]&@@@RandomChoice[iniComplexityZero])/.Rule->List,WorkingPrecision->30]

solsComplexityZero=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->f35,ϵ->ϵϵ},FindRoot[zeroComplexityEqs/.ϵ->ϵϵ,Drop[sol,2]/.Rule->List,WorkingPrecision->20]]],
Join[{f->f35,ϵ->Eini},iniComplexityZero2],Eini-Rest@Range[0,1,(10^(-4))],0<=x<=1&&0<=q1<=1&&Norm[zeroComplexityEqs]<(10^(-6))//.#&]];

Length[solsComplexityZero]

transitionEqs=Simplify[{D[𝒮[f35][ϵ,μ][q1,x],x]/.x->1,SeriesCoefficient[D[𝒮[f35][ϵ,μ][q1,x],q1],{x,1,1}]}];

transitionIni=FindRoot[transitionEqs/.ϵ->Eini,{{μ,μ0[f35][Eini]},{q1,1-0.1RandomReal[]}},WorkingPrecision->20]

solsTransition=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->f35,ϵ->ϵϵ,x->1},FindRoot[transitionEqs/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->f35,ϵ->Eini,x->1},transitionIni],Eini-Range[0,2,(10^(-4))],0<=q1<=1&&𝒮[f35][ϵ,μ][q1,x]>=0/.#&]];

transitionEnd=FindRoot[Append[transitionEqs,𝒮[f35][ϵ,μ][q1,1]],Delete[solsTransition[[-1]],{{1},{3}}]/.Rule->List,WorkingPrecision->50,MaxIterations->500]

(ϵ/.transitionEnd)-E1RSBp[f35][[1]]

zoomPlot=Show[
ListPlot[{{ϵ,(10^3)(μ-μ0[f35][ϵ])}/.solsTransition,{ϵ,(10^3)(μ-μ0[f35][ϵ])}/.solsComplexityZero,{{-1.3,0},{0.1,0}}},Joined->True,PlotRange->{{-1.25,-0.85},{-0.25,1.15}},Epilog->{EdgeForm[Black],Red,Disk[{#,0},Offset[2]]&/@E1RSBp[f35]},PlotStyle->{Red,Black,{Black,Dashed}},Filling->{2->{1}},Prolog->{Red,Opacity[0.3],Polygon[Join[{ϵ,(10^3)(μ-μ0[f35][ϵ])}/.solsTransition,{ϵ,(10^3)(μ-μ0[f35][ϵ])}/.Reverse@solsComplexityZero]]},Frame->True,FrameStyle->Black,LabelStyle->Directive[Black,FontFamily->"Times",FontSize->10],FrameLabel->{None,"(μ - (Subscript[μ, 0])) × (10^3)"},ImageSize->140,FrameTicks->{{{0,1},None},{{-1.2,-1.1,-1,-0.9},None}},Axes->None]
]

nonZoomPlot=ListPlot[{{ϵ,μ}/.solsTransition,{ϵ,μ}/.solsComplexityZero},Joined->True,PlotStyle->Red,Epilog->{Red,Point[{#,μ0[f35][#]}]&/@E1RSBp[f35]}]

rp=RegionPlot[Evaluate[Σa[f35][ϵ,μ]>0&&μ>μm[f35]],{ϵ,-1.5,-1.1},{μ,4.5,5.5},PlotPoints->100,PlotStyle->{Opacity[0.5],Gray},BoundaryStyle->None]

plot35=Show[ContourPlot[0==Σa[f35][ϵ,μ],{ϵ,-1.3,-0.8},{μ,2,5.5},ContourStyle->{{Black,Dashed}}],rp,nonZoomPlot,ImageSize->270,Prolog->Inset[zoomPlot,{-1.13,2.7}],LabelStyle->Directive[Black,FontFamily->"Times",FontSize->10],FrameStyle->Black,FrameLabel->{"E",μ},RotateLabel->False,Exclusions->None,Epilog->{{EdgeForm[Black],Red,Disk[{#,μ0[f35][#]},Offset[2]]&/@E1RSBp[f35]},Line[{{Emax[f]+0.01/.Sqrt[x_]:>-Sqrt[x]/.f->f35,μm[f35]},{-1.08,μm[f35]}}],Text[Style["(Subscript[μ, m])",Directive[Black,FontFamily->"Times",FontSize->10]],{-1.055,μm[f35]},FormatType->"StandardForm"],Line[{{Emax[f35],μm[f35]-0.05},{Emax[f35],4.75}}],Text[Style["(Subscript[E, min])",Directive[Black,FontFamily->"Times",FontSize->10]],{Emax[f35],4.6},FormatType->"StandardForm"],Line[{{E1RSBp[f35][[1]],μ0[f35][E1RSBp[f35][[1]]]-0.05},{E1RSBp[f35][[1]],4.4}}],Text[Style["(Subscript[(SuperMinus[E]),   1RSB])",Directive[Black,FontFamily->"Times",FontSize->10,ScriptMinSize->7]],{E1RSBp[f35][[1]],4.25},FormatType->"StandardForm"],Line[{{E1RSBp[f35][[2]],μ0[f35][E1RSBp[f35][[2]]]-0.05},{E1RSBp[f35][[2]],2.85}}],Text[Style["(Subscript[(SuperPlus[E]),   1RSB])",Directive[Black,FontFamily->"Times",FontSize->10,ScriptMinSize->7],FontVariations->{"CapsType"->"SmallCaps"}],{E1RSBp[f35][[2]],2.7},FormatType->"StandardForm"]}]

Export["figs/complexity_35.pdf",plot35];

Ranges figure

fTest=Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))]/.λ->650/1000;

Eini=E1RSBp[fTest][[2]]-1/10000;

transitionEqs=Simplify[{D[𝒮[fTest][ϵ,μ][q1,x],x]/.x->1,SeriesCoefficient[D[𝒮[fTest][ϵ,μ][q1,x],q1],{x,1,1}]}];

transitionIni=FindRoot[transitionEqs/.ϵ->Eini,{{μ,μ0[fTest][Eini]},{q1,1-0.1RandomReal[]}},WorkingPrecision->20]

solsTransition=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest,ϵ->ϵϵ,x->1},FindRoot[transitionEqs/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest,ϵ->Eini,x->1},transitionIni],Eini-Range[0,2,(10^(-4))],0<=q1<=1&&𝒮[fTest][ϵ,μ][q1,x]>=0/.#&]];

transitionEnd=FindRoot[Append[transitionEqs,𝒮[fTest][ϵ,μ][q1,1]],Delete[solsTransition[[-1]],{{1},{3}}]/.Rule->List,WorkingPrecision->50,MaxIterations->500]

E1RSBp[fTest][[1]]-ϵ/.transitionEnd

fλ=Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))];

transitionEqsλ=Simplify[{D[𝒮[fλ][ϵ,μ][q1,x],x]/.x->1,SeriesCoefficient[D[𝒮[fλ][ϵ,μ][q1,x],q1],{x,1,1}]}];

solsEnd1=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],Prepend[transitionEnd,λ->65/100],Rest[Range[65/100,9/10,(10^(-2))/10]]];

solsEnd2=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],solsEnd1[[-1]],Rest[Range[9/10,945/1000,(10^(-3))/5]]];

solsEnd3=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],solsEnd2[[-1]],Rest[Range[945/1000,960/1000,(10^(-5))]]];

solsEnd4=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],Prepend[transitionEnd,λ->65/100],Rest[Range[65/100,25/100,-(10^(-2))/10]]];

solsEnd5=FoldList[Function[{sol,λλ},Prepend[FindRoot[Append[transitionEqsλ,𝒮[fλ][ϵ,μ][q1,1]]/.λ->λλ,Rest[sol]/.Rule->List,WorkingPrecision->50,MaxIterations->2000],λ->λλ]],solsEnd4[[-1]],Rest[Range[25/100,24/100,-(10^(-4))]]];

solsEndAll=Join[Reverse[solsEnd5],Reverse[solsEnd4],solsEnd1,solsEnd2,solsEnd3];

ListLogPlot[{λ,Abs[ϵ-Emax[fλ]]}/.solsEndAll,PlotRange->All,Joined->True]

Plot[Evaluate[{Ealg3s[s],E1RSBp[Function[q,(1/2)(λ (q^3)+(1-λ)(q^s))]][[1]]}/.s->11],{λ,0,1},PlotPoints->50]

rangePlots=(s|->Show[
Plot[Evaluate[Append[E1RSBp[f],Emax[f]]/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^s))]],{λ,0,1},PlotStyle->{Red,Red,Black},Filling->{1->{2}},PlotRange->{{0,1},{-1.75,1.75}},Frame->True,FrameStyle->Black,LabelStyle->Directive[If[s==5,Black,Opacity[0]]
,FontFamily->"Times",FontSize->10],FrameTicksStyle->{{{},{}},{FontOpacity->0,{}}},FrameLabel->{Style["λ",Opacity[0]],"E"}],
ImageSize->{Automatic,110},ImagePadding->{{45,10},{Automatic,3}},Epilog->Text[Style["s = "<>ToString[s],Directive[Black,FontFamily->"Times",FontSize->10]],{0.8,1}]
])/@{5,8,11,14}

rangePlotsLog=(s|->Show[
LogPlot[Evaluate[Max[#,(10^(-10))]&/@(E1RSBp[f]-Emax[f]/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^s))])],{λ,0,1},PlotStyle->{Red,Red},Filling->{1->{2}},PlotRange->{{0,1},{8 (10^(-8)),8}},Frame->True,FrameStyle->Black,FrameLabel->{λ,If[s==5,"E - (Subscript[E, min])",Style["E - (Subscript[E, min])",Opacity[0]]]},LabelStyle->Directive[Black
,FontFamily->"Times",FontSize->10]],FrameTicksStyle->{{If[s!=5,FontOpacity->0,{}],{}},{{},{}}},ImageSize->{Automatic,110},ImagePadding->{{45,10},{Automatic,3}}
])/@{5,8,11,14}

Table[Export["figs/range_plot_"<>ToString[i]<>".pdf",rangePlots[[i]]],{i,1,3}];
Table[Export["figs/range_plot_log_"<>ToString[i]<>".pdf",rangePlotsLog[[i]]],{i,1,3}];

rangePlotsLogμ=(s|->Show[
LogPlot[Evaluate[Max[#,(10^(-10))]&/@({μm[f]-μ0[f][E1RSBp[f][[2]]],(μm[f]-μ0[f][E1RSBp[f][[1]]])If[μm[f]-μ0[f][E1RSBp[f][[1]]]>2.5,0,1]}/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^s))])],{λ,0,1},PlotStyle->{Red,Red},Filling->{1->{2}},PlotRange->{(10^(-6)),50}]
])/@{5,8,11,14}

lims={FindRoot[Emax[f]-E1RSBp[f][[1]]/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))],{λ,0.83},WorkingPrecision->25],FindRoot[Emax[f]-E1RSBp[f][[1]]/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))],{λ,0.98},WorkingPrecision->25]}

rangeTable14lower=Table[{λ,(E1RSBp[f][[1]]-Emax[f])If[(lims[[1,1,2]])<λ<(lims[[2,1,2]]),0,1]}/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))],{λ,0,1,1/1000}];

rangeTable14upper=Table[{λ,(E1RSBp[f][[2]]-Emax[f])}/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))],{λ,0,1,1/1000}];

solsEndTrimmed=Select[solsEndAll,ϵ<E1RSBp[Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))]][[1]]/.#&];

rangeLowerNew=Join[Select[rangeTable14lower,(λ/.solsEndTrimmed[[1]])>#[[1]]&],{λ,If[0.75<λ<(lims[[2,1,2]]),0,1](ϵ-Emax[Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))]])}/.solsEndTrimmed,Select[rangeTable14lower,(λ/.solsEndTrimmed[[-1]])<#[[1]]&]];

rangeTable14lowerReg=Table[{λ,If[(lims[[1,1,2]])<λ<(lims[[2,1,2]]),Emax[f],E1RSBp[f][[1]]]}/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))],{λ,0,1,1/1000}];

rangeTable14upperReg=Table[{λ,(E1RSBp[f][[2]])}/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))],{λ,0,1,1/1000}];

rangePlot4=Show[
ListPlot[{rangeTable14lowerReg,rangeTable14upperReg},PlotStyle->{Red,Red,{Red,Dashed}},Filling->{1->{2}},PlotRange->{{0,1},{-1.75,1.75}},Frame->True,FrameStyle->Black,LabelStyle->Directive[Opacity[0],FontFamily->"Times",FontSize->10],FrameLabel->{λ,"E"},Joined->True],Plot[Evaluate[Emax[f]/.f->Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))]],{λ,0,1},PlotStyle->Black],ImageSize->{Automatic,110},ImagePadding->{{45,10},{Automatic,3}},Epilog->Text[Style["s = 14",Directive[Black,FontFamily->"Times",FontSize->10]],{0.8,1}]
]

auff2rsb141=FindRoot[{h11Auff[Function[q,(λ (q^3)+(1-λ)(q^14))],x],h12Auff[Function[q,(λ (q^3)+(1-λ)(q^14))],x]},MinimalBy[λ12,Abs[s-14]/.#&][[1,2;;]]/.Rule->List,MaxIterations->1000,WorkingPrecision->20]

auff2rsb142=FindRoot[Ψ[λ,3,14],MinimalBy[λ21,Abs[s-14]/.#&][[1,2;;]]/.Rule->List,MaxIterations->1000,WorkingPrecision->10]

rangePlotLog4=Show[
ListLogPlot[{{#1,Max[If[0.75<#1<0.9,(10^(-10)),#2],(10^(-10))]}&@@@rangeLowerNew,rangeTable14upper,Select[rangeTable14lower,(λ/.solsEndTrimmed[[1]])<#[[1]]<0.9&]},PlotStyle->{Red,Red,{Red,Dashed}},Filling->{1->{2}},PlotRange->{{0,1},{8 (10^(-8)),8}},Frame->True,FrameStyle->Black,LabelStyle->Directive[Black,FontFamily->"Times",FontSize->10],FrameLabel->{λ,Style["E - (Subscript[E, min])",Opacity[0]]},Joined->True,FrameTicksStyle->{{FontOpacity->0,{}},{{},{}}}],ImageSize->{Automatic,110},ImagePadding->{{45,10},{Automatic,3}},Prolog->{{Black,Line[{{λ,-20},{λ,10}}]/.auff2rsb141,Line[{{λ,-20},{λ,10}}]/.auff2rsb142},{Arrowheads[{-0.04,0.04}],Arrow[{{λ/.auff2rsb141,-7},{λ/.auff2rsb142,-7}}]},Text[Style["RSB",Directive[FontFamily->"Times",Black,FontSize->9]],{(((λ/.auff2rsb142)+(λ/.auff2rsb141))/2),-3}],Text[Style["T = 0",Directive[FontFamily->"Times",Black,FontSize->9]],{(((λ/.auff2rsb142)+(λ/.auff2rsb141))/2),-5.5}]}
]

Export["figs/range_plot_4.pdf",rangePlot4];

Export["figs/range_plot_log_4.pdf",rangePlotLog4];

Non-matching limits

fTest1=Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))]/.λ->670/1000;

Eini11=E1RSBp[fTest1][[2]]-1/10000;

transitionEqs1=Simplify[{D[𝒮[fTest1][ϵ,μ][q1,x],x]/.x->1,SeriesCoefficient[D[𝒮[fTest1][ϵ,μ][q1,x],q1],{x,1,1}]}];

transitionIni11=FindRoot[transitionEqs1/.ϵ->Eini11,{{μ,μ0[fTest1][Eini11]},{q1,1-0.1RandomReal[]}},WorkingPrecision->20]

solsTransition11=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest1,ϵ->ϵϵ,x->1},FindRoot[transitionEqs1/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest1,ϵ->Eini11,x->1},transitionIni11],Eini11-Range[0,2,(10^(-3))],0<=q1<=1&&𝒮[fTest1][ϵ,μ][q1,x]>=0/.#&]];

transitionEnd1=FindRoot[Append[transitionEqs1,𝒮[fTest1][ϵ,μ][q1,1]],Delete[solsTransition11[[-1]],{{1},{3}}]/.Rule->List,WorkingPrecision->50,MaxIterations->500]

solsTransition11back=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest1,ϵ->ϵϵ,x->1},FindRoot[transitionEqs1/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest1,x->1},transitionEnd1],(ϵ/.transitionEnd1)+Range[0, (10^(-3)),(10^(-5))],0<=q1<=1&&𝒮[fTest1][ϵ,μ][q1,x]>=0/.#&]];

Eini12=E1RSBp[fTest1][[1]]+(10^(-5));

transitionIni12=FindRoot[transitionEqs1/.ϵ->Eini12,{{μ,μ0[fTest1][Eini12]},{q1,1-0.001RandomReal[]}},WorkingPrecision->20]

solsTransition12=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest1,ϵ->ϵϵ,x->1},FindRoot[transitionEqs1/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest1,ϵ->Eini12,x->1},transitionIni12],Eini12+Range[0, (10^(-3)),(10^(-6))],0<=q1<=1&&𝒮[fTest1][ϵ,μ][q1,x]>=0/.#&]];

zeroComplexity11=Most[FoldWhileList[Function[{sol,μμ},Join[{f->fTest1,μ->μμ},FindRoot[Append[D[𝒮[fTest1][ϵ,μμ][q1,x],{{q1,x}}],𝒮[fTest1][ϵ,μμ][q1,x]],Drop[sol,2]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest1},transitionEnd1[[{2}]],transitionEnd1[[{1,3}]],{x->0.999}],(μ/.transitionEnd1)-Rest@Range[0,5 (10^(-2)),(10^(-4))],0<=q1<=1&&0<=x<=1/.#&]];

zeroComplexity12=Most[FoldWhileList[Function[{sol,μμ},Join[{f->fTest1,μ->μμ},FindRoot[Append[D[𝒮[fTest1][ϵ,μμ][q1,x],{{q1,x}}],𝒮[fTest1][ϵ,μμ][q1,x]],Drop[sol,2]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest1},solsTransition12[[1,{4}]],solsTransition12[[1,{2,5}]],{x->0.999}],(Rationalize[μ/.solsTransition12[[1]],(10^(-20))])-Rest@Range[0, (10^(-1)),(10^(-4))],0<=q1<=1&&0<=x<=1/.#&]];

matchPoint=FindRoot[ Join[transitionEqs1/.q1->q11,transitionEqs1/.q1->q12],Join[solsTransition12[[1,{2,4}]]/.Rule->List,{{q11,q1/.solsTransition11[[-1]]},{q12,q1/.solsTransition12[[1]]}}],WorkingPrecision->30]

inter1RSBtrans=FoldWhileList[Function[{sol,μμ},Join[{f->fTest1,μ->μμ},FindRoot[Join[D[𝒮[fTest1][ϵ,μμ][q11,x1],{{q11,x1}}],D[𝒮[fTest1][ϵ,μμ][q12,x2],{{q12,x2}}],{𝒮[fTest1][ϵ,μμ][q12,x2]-𝒮[fTest1][ϵ,μμ][q11,x1]}],Drop[sol,2]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest1},matchPoint[[{2}]],matchPoint[[{1,3,4}]],{x1->1,x2->1}],(μ/.matchPoint)-Rest[Range[0,(10^(-2)),(10^(-4))]],𝒮[fTest1][ϵ,μ][q11,x1]>=0/.#&];

Length@inter1RSBtrans

orderPlot1=Show[
ListPlot[{Prepend[{ϵ,μ-μ0[fTest1][ϵ]}/.Join[Select[solsTransition12,ϵ<(ϵ/.matchPoint)/.#&],{matchPoint},Select[solsTransition11back,ϵ>(ϵ/.matchPoint)/.#&]],{E1RSBp[fTest1][[1]],0}],{ϵ,μ-μ0[fTest1][ϵ]}/.Join[{transitionEnd1},Select[solsTransition11back,ϵ<(ϵ/.matchPoint)/.#&],{matchPoint},Select[solsTransition12,ϵ>(ϵ/.matchPoint)/.#&]],Prepend[{ϵ,μ-μ0[fTest1][ϵ]}/.zeroComplexity12,{E1RSBp[fTest1][[1]],0}],{ϵ,μ-μ0[fTest1][ϵ]}/.inter1RSBtrans},Joined->True,PlotStyle->{Red,{Red,Dashed},Black,{Red,Dotted}},Filling->{1->{3}},PlotRange->{{-1.40989,-1.40966},{-0.000125,0.00125}},Axes->False],
Plot[0,{ϵ,-1.41,-1.408},PlotStyle->{Black,Dashed}],
Epilog->{{EdgeForm[Black],Red,Disk[{E1RSBp[fTest1][[1]],0},Offset[2]]},Text[Style["λ = 0.67",Directive[Black,FontFamily->"Times",FontSize->10]],Scaled[{0.15,0.85}]],Text[Style["(Subscript[(SuperMinus[E]),   1RSB])",Directive[Black,FontFamily->"Times",FontSize->10,ScriptMinSize->7]],{E1RSBp[fTest1][[1]],0.00045}],Line[{{E1RSBp[fTest1][[1]],0.000045},{E1RSBp[fTest1][[1]],0.00035}}]},Frame->True,FrameStyle->Black,LabelStyle->Directive[FontFamily->"Times",FontSize->10,Black],ImageSize->270,AspectRatio->1/2,FrameLabel->{Style["E",Opacity[0]],"μ - (Subscript[μ, 0])"}
]

fTest2=Function[q,(1/2)(λ (q^3)+(1-λ)(q^14))]/.λ->690/1000;

transitionEqs2=Simplify[{D[𝒮[fTest2][ϵ,μ][q1,x],x]/.x->1,SeriesCoefficient[D[𝒮[fTest2][ϵ,μ][q1,x],q1],{x,1,1}]}];

Eini21=E1RSBp[fTest2][[2]]-1/10000;

transitionIni21=FindRoot[transitionEqs2/.ϵ->Eini21,{{μ,μ0[fTest2][Eini21]},{q1,1-0.1RandomReal[]}},WorkingPrecision->20]

solsTransition21=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest2,ϵ->ϵϵ,x->1},FindRoot[transitionEqs2/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest2,ϵ->Eini21,x->1},transitionIni21],Eini21-Range[0,3,(10^(-3))],0<=q1<=1&&𝒮[fTest][ϵ,μ][q1,x]>=0/.#&]];

transitionEnd2=FindRoot[Append[transitionEqs2,𝒮[fTest2][ϵ,μ][q1,1]],Delete[solsTransition21[[-1]],{{1},{3}}]/.Rule->List,WorkingPrecision->50,MaxIterations->500]

solsTransition212=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest2,ϵ->ϵϵ,x->1},FindRoot[transitionEqs2/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->25]]],
Join[{f->fTest2,x->1},transitionEnd2],(ϵ/.transitionEnd2)+Rest@Range[0,1.2 (10^(-3)),(10^(-6))],0<=q1<=1&&𝒮[fTest2][ϵ,μ][q1,x]>=0/.#&]];

Eini22=E1RSBp[fTest2][[1]]+(10^(-5));

transitionIni22=FindRoot[transitionEqs2/.ϵ->Eini22,{{μ,μ0[fTest2][Eini22]},{q1,1-0.001RandomReal[]}},WorkingPrecision->20]

solsTransition22=Most[FoldWhileList[Function[{sol,ϵϵ},Join[{f->fTest2,ϵ->ϵϵ,x->1},FindRoot[transitionEqs2/.ϵ->ϵϵ,Drop[sol,3]/.Rule->List,WorkingPrecision->25]]],
Join[{f->fTest2,ϵ->Eini22,x->1},transitionIni22],Eini22+Range[0,1.2 (10^(-3)),(10^(-6))],0<=q1<=1&&𝒮[fTest2][ϵ,μ][q1,x]>=0/.#&]];

zeroComplexity21=Most[FoldWhileList[Function[{sol,μμ},Join[{f->fTest2,μ->μμ},FindRoot[Append[D[𝒮[fTest2][ϵ,μμ][q1,x],{{q1,x}}],𝒮[fTest2][ϵ,μμ][q1,x]],Drop[sol,2]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest2},transitionEnd2[[{2}]],transitionEnd2[[{1,3}]],{x->0.999}],(μ/.transitionEnd2)-Rest@Range[0,5 (10^(-2)),(10^(-4))],0<=q1<=1&&0<=x<=1/.#&]];

zeroComplexity22=Most[FoldWhileList[Function[{sol,μμ},Join[{f->fTest2,μ->μμ},FindRoot[Append[D[𝒮[fTest2][ϵ,μμ][q1,x],{{q1,x}}],𝒮[fTest2][ϵ,μμ][q1,x]],Drop[sol,2]/.Rule->List,WorkingPrecision->20]]],
Join[{f->fTest2},solsTransition22[[1,{4}]],solsTransition22[[1,{2,5}]],{x->0.999}],(Rationalize[μ/.solsTransition22[[1]],(10^(-20))])-Rest@Range[0, (10^(-1)),(10^(-4))],0<=q1<=1&&0<=x<=1/.#&]];

orderPlot2=Show[
ListPlot[{{{-1.5,0},{0,0}},Prepend[{ϵ,μ-μ0[fTest2][ϵ]}/.solsTransition22,{E1RSBp[fTest2][[1]],0}],{ϵ,μ-μ0[fTest2][ϵ]}/.solsTransition212,{ϵ,μ-μ0[fTest2][ϵ]}/.zeroComplexity21},Joined->True,PlotStyle->{{Black,Dashed},{Red,Dashed},Red,Black},Filling->{3->{4}},PlotRange->{{-1.401635,-1.401545},{-0.0000125,0.000125}},Axes->False],
Epilog->{{EdgeForm[Black],Red,Disk[{E1RSBp[fTest2][[1]],0},Offset[2]]},Text[Style["λ = 0.69",Directive[Black,FontFamily->"Times",FontSize->10]],Scaled[{0.85,0.85}]],Text[Style["(Subscript[(SuperMinus[E]),   1RSB])",Directive[Black,FontFamily->"Times",FontSize->10,ScriptMinSize->7]],{E1RSBp[fTest2][[1]],0.000045}],Line[{{E1RSBp[fTest2][[1]],0.0000045},{E1RSBp[fTest2][[1]],0.000035}}]},Frame->True,FrameStyle->Black,LabelStyle->Directive[FontFamily->"Times",FontSize->10,Black],ImageSize->270,AspectRatio->1/2,FrameLabel->{"E","μ - (Subscript[μ, 0])"}
]

Export["figs/order_plot_1.pdf",orderPlot1];

Export["figs/order_plot_2.pdf",orderPlot2];