Complexity definitions

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

𝒟[f_][μ_]:=With[{μm=μm[f]},(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=μm[f]},(2/(μm))((μ/(μm))-UnitStep[(μ^2)-((μm)^2)]RealSign[μ] (Sqrt[((μ^2)/((μm)^2))-1]))]

β0[f_][E0_,μ0_]:=-(((E0+μ0)f'[1]+E0 f''[1])/(f[1](f'[1]+f''[1])-((f'[1])^2)))
r00[f_][E0_,μ0_]:=((μ0 f[1]+E0 f'[1])/(f[1](f'[1]+f''[1])-((f'[1])^2)))
d00[f_][E0_,μ0_]:=(1/(f'[1]))-((r00[f][E0,μ0])^2)

𝒮1[f_][E_,μ_]:=-(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))))

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

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

groundState[f_,E_,μ_,Σ0_:0,min_:True]:=FindRoot[{𝒮1[f][E,μ]+𝒟[f][μ]-Σ0,D[𝒮1[f][E,μ],μ]+d𝒟[f][μ]},{{E,If[min,1,-1]Eopt[f,μm[f]]},{μ,If[min,1.1,-1.1]μm[f]}},WorkingPrecision->20]
minμState[f_,E_,μ_,Σ0_:0,min_:True]:=FindRoot[{𝒮1[f][E,μ]+𝒟[f][μ]-Σ0,D[𝒮1[f][E,μ],E]},{{E,If[min,1,-1]Eopt[f,μm[f]]},{μ,If[min,1.1,-1.1]μm[f]}}]

𝒮12[f_][E0_,μ0_,E1_,μ1_,q_][β1_,r11d_,r110_,r01_,r10_,q0_]:=With[{β0=β0[f][E0,μ0],r00=r00[f][E0,μ0]},
-(1/2)+β1 E1-r11d μ1+β1(r11d f'[1]-r110 f'[q0])+β0 β1 f[q]+(β0 r01+β1 r10+r00 r01)f'[q]+((r11d-r110)/(1-q0))(r10-q r00)f'[q]+(1/2)(((β1)^2)(f[1]-f[q0])+(r11d^2)f''[1]+2r01 r10 f''[q]-(r110^2)f''[q0]+((r01^2)-((r11d-r110)/(1-q0))(2q r01-(((1-(q^2))r110-(q0-(q^2))r11d)/(1-q0))))(f'[1]-f'[q0])+((1-(q^2))/(1-q0))+((((r10-q r00))^2)/(1-q0))f'[1]-(1/(f'[1]))((((f'[1])^2)-((f'[q])^2))/(f'[1]-f'[q0]))+((r11d-r110)/(1-q0))(r11d f'[1]-r110 f'[q0])+Log[((1-q0)/(f'[1]-f'[q0]))]
)
]

Eβopt[f_,E0_,μ0_,q_][r11d_,r110_,r10_,q0_]:=With[{β0=β0[f][E0,μ0]},
-β0 f[q]-r11d f'[1]-r10 f'[q]+r110 f'[q0]
]

matB[f_][E0_,μ0_,E1_,μ1_,q_][β1_,r110_,r10_,q0_]:=With[{β0=β0[f][E0,μ0],r00=r00[f][E0,μ0]},
{
{β0 f''[q]+r10 f'''[q],f''[q],0,0,0},
{f''[q],0,0,0,0},
{0,0,-β1 f''[q0]-r110 f'''[q0],-f''[q0],0},
{0,0,-f''[q0],0,0},
{0,0,0,0,0}
}
]

matC[f_][E0_,μ0_,E1_,μ1_,q_][β1_,r11d_,r110_,r01_,r10_,q0_]:=Module[
{Ce,d00=d00[f][E0,μ0],r00=r00[f][E0,μ0]},
Ce[1,1]=d00 f'[1];
Ce[1,2]=r00 f'[1];
Ce[2,2]=-f'[1];
Ce[1,3]=((r11d-r110)/(1-q0))(r01-q ((r11d-r110)/(1-q0)))(f'[1]-f'[q0])+(1/(1-q0))(q f'[1]d00+r00(r10 f'[1]+(r11d-r110)f'[q]));
Ce[1,5]=r00 f'[q]+(r01-q ((r11d-r110)/(1-q0)))(f'[1]-f'[q0]);
Ce[1,4]=-Ce[1,5];
Ce[2,3]=((q r00-r10)/(1-q0))f'[1]-((r11d-r110)/(1-q0))f'[q];
Ce[2,4]=f'[q];
Ce[2,5]=-Ce[2,4];
Ce[3,3]=-((1-(q^2))/(((1-q0))^2))-((r11d-r110)/(1-q0))(((r11d-r110)/(1-q0))f'[1]-2(((q r01-r110)/(1-q0))+((1-(q^2))/(1-q0))((r11d-r110)/(1-q0)))(f'[1]-f'[q0])-2 ((q r00-r10)/(1-q0))f'[q])-((((r10-q r00))^2)/(((1-q0))^2))f'[1];
Ce[3,4]=-(q r01-r110)((f'[1]-f'[q0])/(1-q0))-((r11d-r110)/(1-q0))(((1-(q^2))/(1-q0))(f'[1]-f'[q0])-f'[q0])-f'[q] ((q r00-r10)/(1-q0));
Ce[3,5]=-Ce[3,4]-((r11d-r110)/(1-q0))(f'[1]-f'[q0]);
Ce[4,4]=f'[1]-2f'[q0];
Ce[4,5]=f'[q0];
Ce[5,5]=-f'[1];
Array[Ce,{5,5}]/.Ce[i_,j_]:>Ce[j,i]/;i>j
]

VRS[f_][E0_,μ0_,q_]:=With[
{
β0=β0[f][E0,μ0],
r00=r00[f][E0,μ0]
},
-β0 f[q]-r00 f'[q] q-(Sqrt[(1-(q^2))(1-(((f'[q])^2)/((f'[1])^2)))])
]

V1RSB[f_][E0_,μ0_,q_][y1_,z_,q0_]:=With[
{
β0=β0[f][E0,μ0],
r00=r00[f][E0,μ0]
},
-β0 f[q]-r00 f'[q] q-(1/2) (-((y1((q^2)-q0))/(1+y1 z(1-q0)))+ ((f'[1])/y1)+z (f[1]-f[q0])+(1/z)Log[1+z y1(1-q0)]-(((1+z y1(1-q0) )((f'[q])^2))/(y1 f'[1])))
]

paramSaddles[f_][E0_,μ0_,E1_,μ1_,q_][β1_,r11d_,r110_,r01_,r10_,q0_]:=Evaluate[Simplify[D[𝒮12[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{{β1,r11d,r110,r01,r10,q0}}]]]

μ1Saddle[f_][E0_,μ0_,E1_,μ1_,q_][β1_,r11d_,r110_,r01_,r10_,q0_]:=Evaluate[Simplify[D[𝒮12[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],μ1]+d𝒟[f][μ1]]];

E1Saddle[f_][E0_,μ0_,E1_,μ1_,q_][β1_,r11d_,r110_,r01_,r10_,q0_]:=Evaluate[Simplify[D[𝒮12[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],E1]]];

eigenEquation[f_][E0_,μ0_,E1_,μ1_,q_][β1_,r11d_,r110_,r01_,r10_,q0_][yguess_]:=(y/.FindRoot[Eigenvalues[matB[f][E0,μ0,E1,μ1,q][β1,r110,r10,q0]-y matC[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0]][[2]],{y,yguess}])/;AllTrue[{E0,μ0,E1,μ1,q,β1,r11d,r110,r01,r10,q0,yguess},NumericQ]

eigenvalue[f_][E0_,μ0_,E1_,μ1_,q_][β1_,r11d_,r110_,r01_,r10_,q0_][yguess_]:=With[{y=eigenEquation[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][yguess]},y+(1/y)f''[1]]

f23=Function[x,(1/2)((x^2)+(x^3))];

f34=Function[x,(1/2)((x^3)+(x^4))];

f312=Function[x,(1/2)((829/1000)(x^3)+(1-(829/1000))(x^12))];

Solution for the r variables

rSol[f_][E0_,μ0_,E1_,μ1_,q_][β1_,q0_]=FullSimplify[First[Solve[0==D[𝒮12[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{{r11d,r110,r01,r10}}],{r11d,r110,r01,r10}]]]

𝒮12r[f_][E0_,μ0_,E1_,μ1_,q_][β1_,q0_]=Simplify[𝒮12[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0]/.rSol[f][E0,μ0,E1,μ1,q][β1,q0]];

matCr=Simplify[Simplify[FullSimplify[matC[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0]]/.r11d-r110->FullSimplify[r11d-r110/.rSol[f][E0,μ0,E1,μ1,q][β1,q0]]]/.rSol[f][E0,μ0,E1,μ1,q][β1,q0]];

matBr=Simplify[FullSimplify[matB[f][E0,μ0,E1,μ1,q][β1,r110,r10,q0]]/.rSol[f][E0,μ0,E1,μ1,q][β1,q0]];

Eβoptr=Simplify[Eβopt[f,E0,μ0,q][r11d,r110,r10,q0]/.rSol[f][E0,μ0,E1,μ1,q][0,q0]];

β1Saddle[f_][E0_,μ0_,E1_,μ1_,q_][β1_,q0_]=Simplify[Derivative[1,0][𝒮12r[f][E0,μ0,E1,μ1,q]][β1,q0]];

μ1Saddle[f_][E0_,μ0_,E1_,μ1_,q_][β1_,q0_]=Simplify[D[𝒮12r[f][E0,μ0,E1,μ11,q][β1,q0],μ11]+d𝒟[f][μ1]/.μ11->μ1];

E1Saddle[f_][E0_,μ0_,E1_,μ1_,q_][β1_,q0_]=Simplify[D[𝒮12r[f][E0,μ0,E11,μ1,q][β1,q0],E11]/.E11->E1];

(* Precompute for specific functions f *)
q0Saddle[f34][E0_,μ0_,E1_,μ1_,q_][β1_,q0_]=Simplify[Derivative[0,1][𝒮12r[f34][E0,μ0,E1,μ1,q]][β1,q0]];

Numeric definitions

solLineMin[ff_,E00_,μ00_,Σ0_:0,min_:True,Δq_:(10^(-2))]:=Module[
{
groundProps=groundState[ff,E1,μ1,Σ0,min],
eqs=Join[paramSaddles[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{μ1Saddle[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0],𝒮12[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0]+𝒟[ff][μ1]-Σ0}],
ini
},
ini=Join[{f->ff,E0->E00,μ0->μ00,q->0,β1->β0[ff][E1,μ1]/.groundProps,r11d->r00[ff][E1,μ1]/.groundProps,r110->0,r01->0,r10->0,q0->0},groundProps];
Most[FoldWhileList[
Function[{sol,qq},Check[Join[Take[ini,3],{q->qq},FindRoot[eqs/.q->qq,Drop[sol,4]/.Rule->List,WorkingPrecision->20]],{}]],
ini,
Most[Rest[Range[0,1,Δq]]],
Length[#]>0&
]]]

solLineMinμ[ff_,E00_,μ00_,Σ0_:0,min_:True,Δq_:(10^(-2))]:=Module[
{
groundProps=minμState[ff,E1,μ1,Σ0,min],
eqs=Join[paramSaddles[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{E1Saddle[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0],𝒮12[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0]+𝒟[ff][μ1]-Σ0}],
ini
},
ini=Join[{f->ff,E0->E00,μ0->μ00,q->0,β1->β0[ff][E1,μ1]/.groundProps,r11d->r00[ff][E1,μ1]/.groundProps,r110->0,r01->0,r10->0,q0->0},groundProps];
Most[FoldWhileList[
Function[{sol,qq},Check[Join[Take[ini,3],{q->qq},FindRoot[eqs/.q->qq,Drop[sol,4]/.Rule->List,WorkingPrecision->20]],{}]],
ini,
Most[Rest[Range[0,1,Δq]]],
Length[#]>0&
]]]

solLineMarginal[ff_,E00_,μ00_,Σ0_:0,min_:True,Δq_:(10^(-2))]:=Module[
{
groundProps=FindRoot[𝒮1[ff][E1,μm[ff]]+𝒟[ff][μm[ff]]-Σ0,{E1,Eopt[ff,μm[ff]]-If[min,1,-1]5 (10^(-2))}],
eqs=Simplify@Join[D[𝒮12r[ff][E00,μ00,E1,μm[ff],q][β1,q0],{{β1,q0}}],{𝒮12r[ff][E00,μ00,E1,μm[ff],q][β1,q0]+𝒟[ff][μm[ff]]-Σ0}],
ini
},
ini=Join[{f->ff,E0->E00,μ0->μ00,μ1->μm[ff],q->0,β1->β0[ff][E1,μm[ff]]/.groundProps,q0->0},groundProps];
Most[FoldWhileList[
Function[{sol,qq},Check[Join[Take[ini,4],{q->qq},FindRoot[eqs/.q->qq,Drop[sol,5]/.Rule->List,WorkingPrecision->20]],{}]],
ini,
Most[Rest[Range[0,1,Δq]]],
Length[#]>0&
]]]

solLineConst[ini_,E1f_,ΔE1_]:=Module[
{
eqs=Simplify[Join[D[𝒮12r[f][E0,μ0,E11,μ11,q][β11,q00],{{β11,q00}}],{D[𝒮12r[f][E0,μ0,E11,μ11,q][β11,q00],μ11]+d𝒟[f][μ11]}]/.ini/.{E11->E1,μ11->μ1,β11->β1,q00->q0}]
},
FoldWhileList[
Function[{sol,E11},Check[Join[Take[ini,4],{E1->E11},FindRoot[eqs/.E1->E11,Drop[sol,5]/.Rule->List,WorkingPrecision->20]],{}]],
Insert[DeleteCases[ini,Rule[E1,_]],E1->(E1/.ini),5],
Range[E1/.ini,E1f,ΔE1],
Length[#]>0&
]]

solLineDominantFixedμ1[ff_,E00_,μ00_,μ11_,Δq_:(10^(-2))]:=Module[
{
gs=FindRoot[d𝒟[ff][μ11]+Derivative[0,1][𝒮1[ff]][E1,μ11],{E1,Eopt[ff,μm[ff]]}],
eqs=Join[paramSaddles[ff][E00,μ00,E1,μ11,q][β1,r11d,r110,r01,r10,q0],{μ1Saddle[ff][E00,μ00,E1,μ11,q][β1,r11d,r110,r01,r10,q0]}],
ini
},
ini={f->ff,E0->E00,μ0->μ00,μ1->μ11,q->0,E1->(E1/.gs),β1->β0[ff][E1/.gs,μ11],r11d->r00[ff][E1/.gs,μ11]/.gs,r110->0,r01->0,r10->0,q0->0};
Most@FoldWhileList[
Function[{sol,qq},Check[Join[Take[ini,4],{q->qq},FindRoot[eqs/.q->qq,Drop[sol,5]/.Rule->List,WorkingPrecision->20]],{}]],
ini,
Rest@Most[Range[0,1,Δq]],
Length[#]>0&
]]

eigLineDominantFixedE1[ini0_,δE1_:(10^(-1)),ΔE1_:(10^(-2))]:=Module[
{
eqs=Join[paramSaddles[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{μ1Saddle[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],eigenEquation[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][μm[f]]-μm[f]/2}]/.Take[ini0,3]
},
Rest[FoldWhileList[
Function[{sol,E11},Join[Take[sol,3],FindRoot[eqs/.E1->E11,Most[Drop[sol,3]]/.Rule->List,WorkingPrecision->20],{E1->E11}]],
Permute[ini0,Cycles[{{11,12}}]],
Range[E1/.ini0,E1+δE1/.ini0,ΔE1],
Length[#]>0&
]]]

eigLineDominantFixedE12[ini0_,Δμ1_:(10^(-2))]:=Module[
{
eqs=Join[paramSaddles[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{μ1Saddle[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],eigenvalue[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][μm[f]/2+1]-μ1}]/.Take[ini0,3]
},
Rest[FoldWhileList[
Function[{sol,E11},Check[Join[Take[sol,3],FindRoot[eqs/.E1->E11,Most[Drop[sol,3]]/.Rule->List,WorkingPrecision->20],{E1->E11}],{}]],
Permute[ini0,Cycles[{{11,12}}]],
Range[E1/.ini0,Eopt[f,μm[f]]/.ini0,Δμ1],
Length[#]>0&
]]]

solLineDominantFixedE1[ff_,E00_,μ00_,E11_,Δq_:(10^(-2))]:=Module[
{
gs=FindRoot[Derivative[1,0][𝒮1[ff]][E11,μ1],{μ1,μm[ff]}],
eqs=Join[paramSaddles[ff][E00,μ00,E11,μ1,q][β1,r11d,r110,r01,r10,q0],{E1Saddle[ff][E00,μ00,E11,μ1,q][β1,r11d,r110,r01,r10,q0]}],
ini
},
ini={f->ff,E0->E00,μ0->μ00,E1->E11,q->0,μ1->(μ1/.gs),β1->β0[ff][E11,μ1/.gs],r11d->r00[ff][E11,μ1]/.gs,r110->0,r01->0,r10->0,q0->0};
Most[FoldWhileList[
Function[{sol,qq},Check[Join[Take[ini,4],{q->qq},FindRoot[eqs/.q->qq,Drop[sol,5]/.Rule->List,WorkingPrecision->20]],{}]],
ini,
Rest@Most[Range[0,1,Δq]],
Length[#]>0&
]]]

eigLineDominantFixedμ1[ini0_,δμ1_:(10^(-1)),Δμ1_:(10^(-2))]:=Module[
{
eqs=Join[paramSaddles[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{E1Saddle[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],eigenEquation[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][μm[f]/2+1]-μm[f]/2}]/.Take[ini0,3]
},
Rest[FoldWhileList[
Function[{sol,μ11},Join[Take[sol,3],FindRoot[eqs/.μ1->μ11,Most[Drop[sol,3]]/.Rule->List,WorkingPrecision->20],{μ1->μ11}]],
ini0,
Range[μ1/.ini0,μ1-δμ1/.ini0,-Δμ1],
Length[#]>0&
]]]

eigLineDominantFixedμ12[ini0_,Δμ1_:(10^(-2))]:=Module[
{
eqs=Join[paramSaddles[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{E1Saddle[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],eigenvalue[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][μm[f]/1+1]-μ1}]/.Take[ini0,3]
},
Rest[FoldWhileList[
Function[{sol,μ11},Join[Take[sol,3],FindRoot[eqs/.μ1->μ11,Most[Drop[sol,3]]/.Rule->List,WorkingPrecision->20],{μ1->μ11}]],
ini0,
Range[μ1/.ini0,μm[f]/.ini0,-Δμ1],
Length[#]>0&
]]]

solLineDominant[ff_,E00_,μ00_,Δq_:(10^(-2))]:=Module[
{
eqs=Join[paramSaddles[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{E1Saddle[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0],μ1Saddle[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0]}],
ini
},
ini={f->ff,E0->E00,μ0->μ00,q->0,E1->0,μ1->0,β1->β0[ff][0,0],r11d->r00[ff][0,0],r110->0,r01->0,r10->0,q0->0};
FoldWhileList[
Function[{sol,qq},Check[Join[Take[ini,3],{q->qq},FindRoot[eqs/.q->qq,Drop[sol,4]/.Rule->List,WorkingPrecision->20]],{}]],
ini,
Most[Rest[Range[0,1,Δq]]],
Length[#]>0&
]]

franzParisi1RSB[ff_,E00_,μ00_,Δq_:(10^(-3))]:=Module[{
ini=Join[{f->ff,E0->E00,μ0->μ00,q->0,q0->0},FindMaximum[V1RSB[ff][E00,μ00,0][y1,z,0],{{y1,1},{z,1}}][[2]]]
},
FoldList[Function[{sol,qq},Join[{f->ff,E0->E00,μ0->μ00,q->qq},FindRoot[D[V1RSB[ff][E00,μ00,qq][y1,z,q0],{{y1,z,q0}}],Drop[sol,4]/.Rule->List,WorkingPrecision->15]]],
ini,
Range[0,1,Δq]
]]

addEigenvalueInformation[data_]:=Rest@FoldList[Function[{sol,dat},Join[dat,FindRoot[Eigenvalues[Evaluate[Evaluate[matBr-y matCr]/.dat]][[2]],{y,y/.sol}]]],Append[First[data],y->μm[f/.First[data]]+1],Rest[data]]

Data generation

blendList[color1_,color2_][data_,center_,shift_:0]:=Module[{
ci=Position[data,center][[1,1]],
n1,n2,
n
},
n1=ci+shift;
n2=Length[data]-ci-shift;
n=Max[n1,n2];
Join[Table[Blend[{White,color1},x],{x,(n1-1)/n+1/n/2,1/n/2,-1/n}],Table[Blend[{White,color2},x],{x,1/n/2,1/n/2+(n2-1)/n,1/n}]]
]

testMin=solLineMin[f34,Eopt[f34,61/10],61/10,0,True,(10^(-4))];
testMax=solLineMin[f34,Eopt[f34,61/10],61/10,0,False,(10^(-4))];

testμMax=solLineMinμ[f34,Eopt[f34,61/10],61/10,0,True,(10^(-4))];
testμmin=solLineMinμ[f34,Eopt[f34,61/10],61/10,0,False,(10^(-4))];

minEigμ=First@MinimalBy[Rest[testμMax],(#+(1/(#))f''[1]&)@eigenEquation[f34][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][4]/.#&]

eigμ1=eigLineDominantFixedμ1[minEigμ,2 (10^(-1)),(10^(-3))];

minEigμ2=First@MinimalBy[Rest[testμMax],(Abs[#+(1/(#))f''[1]-μ1]&)@eigenEquation[f34][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][4]/.#&]

eigμ2=eigLineDominantFixedμ12[minEigμ2,(10^(-4))];

minEigE=First@MinimalBy[Rest[testMin],(#+(1/(#))f''[1]&)@eigenEquation[f34][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][4]/.#&]

eigE1=eigLineDominantFixedE1[minEigE,3 (10^(-2)),(10^(-3))];

minEigE2=First@MinimalBy[Rest[testMin],(Abs[#+(1/(#))f''[1]-μ1]&)@eigenEquation[f34][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][4]/.#&]

eigE2=eigLineDominantFixedE12[minEigE2,(10^(-4))];

eigμshape1=Polygon[{q,μ1}/.Join[eigμ2,Reverse[Select[eigμ1,(μ1/.#)>(μ1/.eigμ2[[-2]])&]],Select[testμMax,(q/.eigμ1[[1]])<=(q/.#)<=(q/.eigμ2[[1]])&]]];

eigμshape2=Polygon[{q,μ1}/.Join[Reverse@Select[eigμ1,(μ1/.#)<(μ1/.eigμ2[[-2]])&],Reverse@eigμ2,Select[testμMax,(q/.eigμ2[[1]])<=(q/.#)&],Reverse@Select[testμmin,(q/.eigμ2[[1]])<=(q/.#)&]]]

eigEshape1=Polygon[{q,E1}/.Join[Most@eigE2,Reverse[Select[eigE1,(E1/.#)<(E1/.eigE2[[-2]])&]],Select[testMin,(q/.eigE1[[1]])<=(q/.#)<=(q/.eigE2[[1]])&]]]

eigEshape2=Polygon[{q,E1}/.Join[Reverse@Select[eigE1,(E1/.#)<(E1/.eigE2[[-2]])&],Reverse@Most@eigE2,Select[testMin,(q/.eigE2[[1]])<=(q/.#)&],Reverse@Select[testMax,(q/.eigE2[[1]])<=(q/.#)&]]];

E1levels=Range[Eopt[f34,μm[f34]]-3/100,Eopt[f34,μm[f34]]+3/100,1/200];

testElines=solLineDominantFixedE1[f34,Eopt[f34,61/10],61/10,#,(10^(-3))]&/@E1levels;

adjustedElines=((sl|->Flatten[Table[Join[sl[[i]],Select[If[(μ1/.sl[[i,-1]])>(μ1/.Last[Most[testμMax]]),testμMax,testμmin],(q/.sl[[i,-1]])<(q/.#)<(q/.sl[[i+1,1]])&]],{i,1,Length[sl]-1}],1])@Append[SplitBy[#,If[q==0/.#,1,Sign[𝒮12[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0]+𝒟[f][μ1]]/.#]&][[;;;;2]],{{q->1}}])&/@testElines;

plotlegendμ=BarLegend[{x|->styletable[[Ceiling[(Length[styletable]-1)((x-First[E1levels]+1/200)/(E1levels[[-1]]-First[E1levels]))]]],{First[E1levels]-1/200,E1levels[[-1]]}},LabelStyle->Directive[FontFamily->Times,Black,FontSize->9],Ticks->{-1.72,-1.71,-1.70,{Eopt[f34,μm[f34]], Text[Style["(Subscript[E, th])",Directive[FontFamily->Times,Black]],FormatType->"StandardForm"]},-1.68,-1.67},LegendMarkerSize->120,LegendLabel->"(Subscript[E, 1])",FrameStyle->Black]

styletable=blendList[ColorData[97][2],ColorData[97][3]][E1levels,Eopt[f34,μm[f34]]];
p2=Show[
ListPlot[{q,μ1}/.Join[{testμmin,testμMax},adjustedElines],Joined->True,Filling->{Table[i->{{i+1},styletable[[i-1]]},{i,2,Length[testElines]+1}]},PlotRange->{{0,1},{5.93,6.17}},PlotStyle->Opacity[0],LabelStyle->Directive[FontFamily->Times,Black],FrameLabel->{q,"(Subscript[μ, 1])"},FrameStyle->Black,Frame->True],
ListPlot[{q,μ1}/.adjustedElines,PlotStyle->{{Black,Thin}},Joined->True],
ListPlot[{{q,μ1}/.Join[testμmin,Reverse@testμMax],{q,μ1}/.testMin},PlotStyle->{Black,{Black,Thin,DotDashed}},Joined->True],
Epilog->{{Dashed,Line[{{0,μm[f34]},{0.85,μm[f34]}}]},{Directive[FontFamily->Times],Text[Style["(Subscript[μ, m])",Directive[FontFamily->Times,Black]],{0.89,μm[f34]},FormatType->"StandardForm"],Text["(Subscript[μ, 0])",{0.34,61/10}]},{Dotted,Line[{{0,61/10},{0.3,61/10}}]},{Opacity[0.3],eigμshape2},{Opacity[0.2],eigμshape1}},ImageSize->250,RotateLabel->False
]

SetDirectory[NotebookDirectory[]];Export["figs/gapped_min_stability.pdf",p2];Export["figs/gapped_min_stability_legend.pdf",plotlegendμ];

testμMax[[-1]]

testMax[[-1]]

μ1levels=Range[5900/1000,6125/1000,25/1000/2];

testμlines=solLineDominantFixedμ1[f34,Eopt[f34,61/10],61/10,#,(10^(-3))]&/@μ1levels;

adjustedμlines=((sl|->Flatten[Table[Join[sl[[i]],Select[If[(E1/.sl[[i,-1]])>(E1/.Last[Most[testMax]]),testMax,testMin],(q/.sl[[i,-1]])<(q/.#)<(q/.sl[[i+1,1]])&]],{i,1,Length[sl]-1}],1])@Append[SplitBy[#,If[q==0/.#,1,Sign[𝒮12[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0]+𝒟[f][μ1]]/.#]&][[;;;;2]],{{q->1}}])&/@testμlines;

plotlegend=BarLegend[{x|->styletable[[Ceiling[(Length[styletable]-1)((x-First[μ1levels])/(Last[μ1levels]-First[μ1levels]))]]],{First[μ1levels],Last[μ1levels]}},LabelStyle->Directive[FontFamily->Times,Black,FontSize->9],Ticks->{5.90,5.95,{6, Text[Style["(Subscript[μ, m])",Directive[FontFamily->Times,Black]],FormatType->"StandardForm"]},6.05,6.1},LegendMarkerSize->120,LegendLabel->μ,FrameStyle->Black]

styletable=blendList[ColorData[97][1],ColorData[97][4]][Range[5900/1000,6125/1000,25/1000/2],6,-1];
p1=Show[ListPlot[{q,E1}/.Join[{testMin,testMax},adjustedμlines],Joined->True,Filling->{Append[Table[i->{{i+1},styletable[[i-2]]},{i,3,Length[testμlines]+1}],{1->{{Length[testμlines]+2},styletable[[Length[testμlines]]]}}]},PlotRange->{{0,1},{-1.735,-1.675}},PlotStyle->Opacity[0],LabelStyle->Directive[FontFamily->Times,Black],FrameLabel->{q,"(Subscript[E, 1])"},FrameStyle->Black,Frame->True],
ListPlot[{q,E1}/.adjustedμlines,PlotStyle->{{Black,Thin}},Joined->True],
ListPlot[{{q,E1}/.Join[testMin,Reverse@testMax],{q,E1}/.testμMax},PlotStyle->{Black,{Black,Thin,DotDashed}},Joined->True],
Epilog->{{Dashed,Line[{{0,Eopt[f34,μm[f34]]},{0.85,Eopt[f34,μm[f34]]}}]},{Directive[FontFamily->Times],Text["(Subscript[E, th])",{0.89,Eopt[f34,μm[f34]]}],Text["(Subscript[E, 0])",{0.34,Eopt[f34,61/10]}]},{Dotted,Line[{{0,Eopt[f34,61/10]},{0.3,Eopt[f34,61/10]}}]},{Opacity[0.3],eigEshape2},{Opacity[0.2],eigEshape1}},ImageSize->250,RotateLabel->False
]

SetDirectory[NotebookDirectory[]];Export["figs/gapped_min_energy.pdf",p1];Export["figs/gapped_min_energy_legend.pdf",plotlegend];

fptest=franzParisi1RSB[f34,Eopt[f34,61/10],61/10,(10^(-4))];

t1=Interpolation[{q,V1RSB[f][E0,μ0,q][y1,z,q0]}/.Rest@fptest];
t2=Interpolation[{q,E1}/.testMin];

FindMinimum[t2[x]-t1[x],{x,0.6}]

Plot[t2[x]-t1[x],{x,0,0.8},PlotRange->All]

p3=Show[ListPlot[{{q,V1RSB[f][E0,μ0,q][y1,z,q0]}/.fptest,{q,E1}/.testMin},Joined->True,Frame->True,FrameStyle->Black,FrameLabel->{q,"(Subscript[E, 1])"},LabelStyle->Directive[Black,FontFamily->Times,FontSize->9],PlotRange->{Automatic,-1.685}],RotateLabel->False,Epilog->{{Line[{{0.78,Eopt[f34,μm[f34]]},{0.9,Eopt[f34,μm[f34]]}}]},{Directive[FontFamily->Times,FontSize->9],Text["(Subscript[E, th])",{0.95,Eopt[f34,μm[f34]]}],Text["(Subscript[E, 0])",{0.84,Eopt[f34,61/10]}]},{Line[{{0.98,Eopt[f34,61/10]},{0.88,Eopt[f34,61/10]}}]}},ImageSize->250]

SetDirectory[NotebookDirectory[]];Export["figs/franz_parisi.pdf",p3];

Plotting

testMarge=solLineMarginal[f34,Eopt[f34,μm[f34]]+2/1000,μm[f34],0,True,(10^(-3))];
testMargeMax=solLineMarginal[f34,Eopt[f34,μm[f34]]+2/1000,μm[f34],0,False,(10^(-3))];

testMarge2=solLineMarginal[f34,Eopt[f34,μm[f34]]+1/1000,μm[f34],0,True,(10^(-3))];
testMargeMax2=solLineMarginal[f34,Eopt[f34,μm[f34]]+1/1000,μm[f34],0,False,(10^(-3))];

testMarge3=solLineMarginal[f34,Eopt[f34,μm[f34]]+1/5000,μm[f34],0,True,(10^(-4))];
testMargeMax3=solLineMarginal[f34,Eopt[f34,μm[f34]]+1/5000,μm[f34],0,False,(10^(-4))];

testMarge=solLineMarginal[f34,Eopt[f34,μm[f34]]+2/1000,μm[f34],0,True,(10^(-3))];

testDom=solLineDominant[f34,Eopt[f34,μm[f34]]+1/10000,μm[f34],(10^(-3))];

ListLogPlot[{{q,y+(1/y)f''[1]-μm[f]}/.test2,{q,y+(1/y)f''[1]-μ1}/.test2},Joined->True]

test[[1]]

ListPlot[{{q,μ1-μm[f]}/.testMin,{q,μ1-μm[f]}/.testMax,{q,μ1-μm[f]}/.testDom},Joined->True,PlotRange->{-0.01,0.01}]

Show[
ListPlot[{q,V1RSB[f][E0,μ0,q][y1,z,q0]-E0}/.fptest,Joined->True],
Plot[VRS[f][E0,μ0,q]-E0/.Take[fptest[[1]],3],{q,0,1},PlotStyle->Dashed,PlotRange->All]
]

ListLogPlot[{q,Eβoptr-E0}/.testDom,Joined->True]

ListLogPlot[{{q,E1-E0}/.testMarge2,{q,E1-E0}/.testMargeMax2,{q,E1-E0}/.testMarge,{q,E1-E0}/.testMargeMax,{q,E1-E0}/.testMarge3,{q,E1-E0}/.testMargeMax3 ,{q,E1-E0}/.testMarge4,{q,E1-E0}/.testMargeMax4,{q,Eβoptr-E0}/.testDom },Joined->True]

ListLogPlot[{{q,E1-E0}/.testMin,{q,E1-E0}/.testMax,{q,E1-E0}/.testMarge,{q,E1-E0}/.testMargeMax,{q,V1RSB[f][E0,μ0,q][y1,z,q0]-E0}/.fptest },Joined->True]

ListPlot[{q,E1}/.testMin,Joined->True]

ListLogPlot[{{q,μ1-μm[f34]}/.testMin},Joined->True]

ListLogPlot[Differences[μ1/.testMin],Joined->True]