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[Det[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]],{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))];

f39:=q|->(1/2)((781/1000)(q^3)+(1-(781/1000))(q^9))

f38:=Function[x,(1/2)((760/1000)(x^3)+(1-(760/1000))(x^8))]

q12SF[f_]:=(Sqrt[1-((f'[1])/(f''[1]))])

χmg[f_]:=1/(Sqrt[f''[1]])

βSF[f_]:=((q12SF[f])/(χmg[f] f'[q12SF[f]]))

ESF[f_]:=-χmg[f]f'[1]-βSF[f]f[q12SF[f]]

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)),constSteps_:True,minapproach_:4]:=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,
If[constSteps,Most[Rest[Range[0,1,Δq]]],1-10^(-Rest@Range[0,minapproach,Δq])],
Length[#]>0&
]]]

solLineMinμ[ff_,E00_,μ00_,Σ0_:0,min_:True,Δq_:(10^(-2)),constSteps_:True,minapproach_:4]:=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,MaxIterations->1000]],{}]],
ini,
If[constSteps,Most[Rest[Range[0,1,Δq]]],1-10^(-Rest@Range[0,minapproach,Δq])],
Length[#]>0&
]]]

solLineMarginalMin[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]}],
eqs=Join[paramSaddles[ff][E00,μ00,E1,μm[ff],q][β1,r11d,r110,r01,r10,q0],{𝒮12[ff][E00,μ00,E1,μm[ff],q][β1,r11d,r110,r01,r10,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,r11d->r00[ff][E1,μm[ff]]/.groundProps,r110->0,r01->0,r10->0,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&
]]]

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&
]]

solLineDominantFixedμ[ff_,E00_,μ00_,μ11_,Δq_:(10^(-2))]:=Module[
{
eqs=Join[paramSaddles[ff][E00,μ00,E1,μ11,q][β1,r11d,r110,r01,r10,q0],{E1Saddle[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->Eopt[ff,μ11],β1->β0[ff][Eopt[ff,μ11],μ11],r11d->r00[ff][Eopt[ff,μ11],μ11],r110->0,r01->0,r10->0,q0->0};
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&
]]

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->30]]],
ini,
Range[0,1,Δq]
]]

addEigenvalueInformation[data_]:=Map[Append[#,y->(eigenEquation[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][μm[f]/1+1]/.#)]&,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}]]
]

mostCommonData[ff_,E00_,μ00_]:=Module[
{
minEnergyLine=solLineMin[ff,E00,μ00,0,True,(10^(-4))],
maxEnergyLine=solLineMin[ff,E00,μ00,0,False,(10^(-4))],
maxStabilityLine=solLineMinμ[ff,E00,μ00,0,True,(10^(-4))],
minStabilityLine=solLineMinμ[ff,E00,μ00,0,False,(10^(-4))],
eigenIsoStartEnergy,eigenIsoStartStability,eigenMinStartEnergy,eigenMinStartStability,
minE=E1/.minμState[ff,E1,μ1],Eth,Erange,maxμ=μ1/.groundState[ff,E1,μ1],μrange
},
Eth=Eopt[ff,μm[ff]];
Erange=Rationalize[Rest[Range[minE,minE+2(Eth-minE),(Eth-minE)/10]],(10^(-10))];
μrange=Rationalize[Reverse[Rest[Range[maxμ,maxμ-2(maxμ-μm[ff]),-(maxμ-μm[ff])/10]]],(10^(-10))];
eigenIsoStartStability=First@MinimalBy[Rest[maxStabilityLine],eigenvalue[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0][1+μm[ff]]/.#&];
eigenIsoStartEnergy=First@MinimalBy[Rest[minEnergyLine],eigenvalue[ff][E00,μ00,E1,μ1,q][β1,r11d,r110,r01,r10,q0][1+μm[ff]]/.#&];
eigenMinStartStability=First@MinimalBy[Rest[maxStabilityLine],(If[#>μm[ff]/2,Abs[#+(1/(#))f''[1]-μ1],(10^10)]&)@eigenEquation[f34][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][4]/.#&];
eigenMinStartEnergy=First@MinimalBy[Rest[minEnergyLine],(If[#>μm[ff]/2,Abs[#+(1/(#))f''[1]-μ1],(10^10)]&)@eigenEquation[f34][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0][4]/.#&];
<|
"minEnergyLine"->minEnergyLine,
"maxEnergyLine"->maxEnergyLine,
"maxStabilityLine"->maxStabilityLine,
"minStabilityLine"->minStabilityLine,
"eigenIsoStability"->eigLineDominantFixedμ1[eigenIsoStartStability,2 (10^(-1)),(10^(-3))],
"eigenIsoEnergy"->eigLineDominantFixedE1[eigenIsoStartEnergy,3 (10^(-2)),(10^(-3))],
"eigenMinStability"->eigLineDominantFixedμ12[eigenMinStartStability,(10^(-4))],
"eigenMinEnergy"->eigLineDominantFixedE12[eigenMinStartEnergy,(10^(-5))],
"energyIsolines"->(solLineDominantFixedE1[ff,E00,μ00,#,(10^(-3))]&/@Erange),
"stabilityIsolines"->(solLineDominantFixedμ1[ff,E00,μ00,#,(10^(-3))]&/@μrange)
|>
]

adjustEnergyLines[data_]:=((sl|->Flatten[Table[Join[sl[[i]],Select[If[(μ1/.sl[[i,-1]])>(μ1/.Last[Most[data["maxStabilityLine"]]]),data["maxStabilityLine"],data["minStabilityLine"]],(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}}])&/@data["energyIsolines"]

Options[mostCommonPlots]={lowEnergyColor->ColorData[97][2],highEnergyColor->ColorData[97][3],lowStabilityColor->ColorData[97][1],highStabilityColor->ColorData[97][4]};
mostCommonPlots[data_,OptionsPattern[]]:=Module[
{
eigenIsoStabilityShape=Polygon[{q,μ1}/.Join[data["eigenMinStability"],Reverse[Select[data["eigenIsoStability"],(μ1/.#)>(μ1/.data["eigenMinStability"][[-2]])&]],Select[data["maxStabilityLine"],(q/.data["eigenIsoStability"][[1]])<=(q/.#)<=(q/.data["eigenMinStability"][[1]])&]]],
eigenMinStabilityShape=Polygon[{q,μ1}/.Join[Reverse@Select[data["eigenIsoStability"],(μ1/.#)<(μ1/.data["eigenMinStability"][[-2]])&],Reverse@data["eigenMinStability"],Select[data["maxStabilityLine"],(q/.data["eigenMinStability"][[1]])<=(q/.#)&],Reverse@Select[data["minStabilityLine"],(q/.data["eigenMinStability"][[1]])<=(q/.#)&]]],
eigenIsoEnergyShape=Polygon[{q,E1}/.Join[Most@data["eigenMinEnergy"],Reverse[Select[data["eigenIsoEnergy"],(E1/.#)<(E1/.data["eigenMinEnergy"][[-2]])&]],Select[data["minEnergyLine"],(q/.data["eigenIsoEnergy"][[1]])<=(q/.#)<=(q/.data["eigenMinEnergy"][[1]])&]]],
eigenMinEnergyShape=Polygon[{q,E1}/.Join[Reverse@Select[data["eigenIsoEnergy"],(E1/.#)<(E1/.data["eigenMinEnergy"][[-2]])&],Reverse@Most@data["eigenMinEnergy"],Select[data["minEnergyLine"],(q/.data["eigenMinEnergy"][[1]])<=(q/.#)&],Reverse@Select[data["maxEnergyLine"],(q/.data["eigenMinEnergy"][[1]])<=(q/.#)&]]],
adjustedEnergyLines=adjustEnergyLines[data],
adjustedStabilityLines=((sl|->Flatten[Table[Join[sl[[i]],Select[If[(E1/.sl[[i,-1]])>(E1/.Last[Most[data["maxEnergyLine"]]]),data["maxEnergyLine"],data["minEnergyLine"]],(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}}])&/@data["stabilityIsolines"],
Erange=E1/.data["energyIsolines"][[All,1]],μrange=μ1/.data["stabilityIsolines"][[All,1]],Eth=Eopt[f,μm[f]]/.data["minEnergyLine"][[1]],plotlegendμ,plotlegendE,styletableμ,styletableE,μmar=μm[f]/.data["minEnergyLine"][[1]],E00=E0/.data["minEnergyLine"][[1]],Emin=E1/.data["minEnergyLine"][[1]],μmax=μ1/.data["maxStabilityLine"][[1]]
},
styletableE=blendList[OptionValue[lowEnergyColor],OptionValue[highEnergyColor]][Erange,Eth];
styletableμ=blendList[OptionValue[lowStabilityColor],OptionValue[highStabilityColor]][μrange,μmar,-1];
plotlegendμ=BarLegend[{x|->styletableE[[Ceiling[(Length[styletableE]-1)((x-First[Erange])/(Erange[[-1]]-First[Erange]))]]],{First[Erange],Erange[[-1]]}},LabelStyle->Directive[FontFamily->Times,Black,FontSize->9],Ticks->{{Eopt[f34,μm[f34]], Text[Style["\!\(\*SubscriptBox[StyleBox[\"E\",FontSlant->\"Italic\"], \"th\"]\)",Directive[FontFamily->Times,Black]],FormatType->"StandardForm"]}},LegendMarkerSize->120,LegendLabel->"\!\(\*SubscriptBox[StyleBox[\"E\",FontSlant->\"Italic\"], \"1\"]\)",FrameStyle->Black];
plotlegendE=BarLegend[{x|->styletableμ[[Ceiling[(Length[styletableμ]-1)((x-First[μrange])/(Last[μrange]-First[μrange]))]]],{First[μrange],Last[μrange]}},LabelStyle->Directive[FontFamily->Times,Black,FontSize->9],Ticks->{{μmar, Text[Style["(Subscript[μ, m])",Directive[FontFamily->Times,Black]],FormatType->"StandardForm"]}},LegendMarkerSize->120,LegendLabel->μ,FrameStyle->Black];
{
Show[
ListPlot[{q,μ1}/.Join[{data["minStabilityLine"],data["maxStabilityLine"]},adjustedEnergyLines],Joined->True,Filling->{Table[i->{{i+1},styletableE[[i-1]]},{i,2,Length[adjustedEnergyLines]+1}]},PlotRange->{{0,1},{μmar-(1/3)(μmax-μmar),μmax+(1/20)(μmax-μmar)}},PlotStyle->Opacity[0],LabelStyle->Directive[FontFamily->Times,Black],FrameLabel->{q,"(Subscript[μ, 1])"},FrameStyle->Black,Frame->True],
ListPlot[{q,μ1}/.adjustedEnergyLines,PlotStyle->{{Black,Thin}},Joined->True,PlotRange->All],
ListPlot[{{q,μ1}/.Join[data["minStabilityLine"],Reverse@data["maxStabilityLine"]],{q,μ1}/.Join[data["minEnergyLine"],Reverse@data["maxEnergyLine"]]},PlotStyle->{Black,{Black,Thin,DotDashed}},Joined->True],
Epilog->{{Dashed,Line[{{0,μmar},{0.85,μmar}}]},{Directive[FontFamily->Times],Text[Style["(Subscript[μ, m])",Directive[FontFamily->Times,Black]],{0.89,μmar},FormatType->"StandardForm"],Text["\!\(\*SubscriptBox[StyleBox[\"μ\",FontSlant->\"Italic\"], \"0\"]\)",{0.34,61/10}]},{Dotted,Line[{{0,61/10},{0.3,61/10}}]},{Opacity[0.3],eigenIsoStabilityShape},{Opacity[0.2],eigenMinStabilityShape}},ImageSize->250,RotateLabel->False,AspectRatio->1
],
plotlegendμ,
Show[ListPlot[{q,E1}/.Join[{data["minEnergyLine"],data["maxEnergyLine"]},adjustedStabilityLines],Joined->True,Filling->{Append[Table[i->{{i+1},styletableμ[[i-2]]},{i,3,Length[adjustedStabilityLines]+1}],{1->{{Length[adjustedStabilityLines]+2},styletableμ[[Length[adjustedStabilityLines]]]}}]},PlotRange->{{0,1},{ Emin-(1/20)(Eth-Emin),Eth+(1/3)(Eth-Emin)}},PlotStyle->Opacity[0],LabelStyle->Directive[FontFamily->Times,Black],FrameLabel->{q,"(Subscript[E, 1])"},FrameStyle->Black,Frame->True],
ListPlot[{q,E1}/.adjustedStabilityLines,PlotStyle->{{Black,Thin}},Joined->True,PlotRange->All],
ListPlot[{{q,E1}/.Join[data["minEnergyLine"],Reverse@data["maxEnergyLine"]],{q,E1}/.Join[data["minStabilityLine"],Reverse@data["maxStabilityLine"]]},PlotStyle->{Black,{Black,Thin,DotDashed}},Joined->True],
Epilog->{{Dashed,Line[{{0,Eth},{0.85,Eth}}]},{Directive[FontFamily->Times],Text["(Subscript[E, th])",{0.89,Eth}],Text["(Subscript[E, 0])",{0.34,E00}]},{Dotted,Line[{{0,E00},{0.3,E00}}]},{Opacity[0.3],eigenIsoEnergyShape},{Opacity[0.2],eigenMinEnergyShape}},ImageSize->250,RotateLabel->False,AspectRatio->1
],
plotlegendE
}
]

(E1/.groundState[f34,E1,μ1])

(-1.675-Eopt[f34,μm[f34]])/((E1/.groundState[f34,E1,μ1])-Eopt[f34,μm[f34]])

testData=mostCommonData[f34,Eopt[f34,61/10],61/10];

ClearAll[mostCommonPlots]

mostCommonPlots[testData,highStabilityColor->Red]

testData38=mostCommonData[f38,Eopt[f38,61/10],61/10];

ListPlot[{q,μ1}/.testData38["energyIsolines"],Joined->True]

ListPlot[{q,μ1}/.adjustEnergyLines[testData38]]

mostCommonPlots[testData38]

dtest=Evaluate[D[𝒮1[f38][E1,μ1],μ1]]+d𝒟[f38][μ1];

ContourPlot[{0==𝒮1[f38][E1,μ1]+𝒟[f38][μ1],0==dtest},{E1,-1.3,-0.95},{μ1,5.9,6.4},Exclusions->None,PlotPoints->30]

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))];

testDominant=solLineDominant[f34,Eopt[f34,61/10],61/10,(10^(-3))];

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}/.Join[testMin,Reverse@testMax]},PlotStyle->{Black,{Black,Thin,DotDashed},{Black,Thin}},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,AspectRatio->1
]

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}/.Join[testμmin,Reverse@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,AspectRatio->1
]

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

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

ListPlot[{q,V1RSB[f][E0,μ0,q][y1,z,q0]}/.fptest]

𝒱∞2=(1/2) (((((q^2)-q0) y1)/(-1+(-1+q1+q0 x1-q1 x1) y1 z))+z f[1]+2 β0 f[q]-x1 z f[q0]+(-1+x1) z f[q1]+((x1 Log[1+y1 z-q1 y1 z]+Log[((-1+(-1+q1+q0 x1-q1 x1) y1 z)/(-1+(-1+q1) y1 z))])/(x1 z))+((f'[1])/y1)+f'[q] (2 q r0+(((-1+(-1+q1+q0 x1-q1 x1) y1 z) f'[q])/(y1 f'[1]))))/.{r0->r00[f][E0,μ0],β0->β0[f][E0,μ0]};

fpHess=Simplify[D[𝒱∞2,{{y1,x1,q0,q1,z},2}]/.q1->q0]/.x1->1/2//Simplify;

Simplify[D[𝒱∞2,{{y1,x1,q0,q1,z},1}]/.q1->q0]/.x1->1/2/.fptest[[1]]

fptest[[5911]]

𝒱∞2/.

Select[
Table[FindMinimum[𝒱∞2/.fptest[[-5600,;;4]],{{y1,4 RandomReal[]},{x1,RandomReal[],0,1},{q0,RandomReal[],0,1},{q1,RandomReal[],0,1},{z,2RandomReal[]}}],{10000}],
(10^(-10))>Norm[D[𝒱∞2/.fptest[[-5600,;;4]],{{y1,x1,q0,q1,z}}]]&&Abs[q0-q1]>(10^(-5))&&Abs[q1-1]>(10^(-5))/.#[[2]]&
]

dat=Transpose[Eigenvalues[fpHess/.#]&/@fptest];

ListPlot[RealAbs@Chop@dat,Joined->True]

MaximalBy[fptest,V1RSB[f][E0,μ0,q][y1,z,q0]/.#&]

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

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

testMinMarge2=solLineMarginalMin[f38,Eopt[f38,μm[f38]]+#,μm[f38],0,True,(10^(-4))]&/@Range[0,25/1000,5/1000];

testMaxMarge2=solLineMarginalMin[f38,Eopt[f38,μm[f38]]+#,μm[f38],0,False,(10^(-4))]&/@Range[0,25/1000,5/1000];

Eopt[f38,μm[f38]]+Range[-2/100,25/1000,5/1000]//N

testMinMarge=solLineMin[f38,Eopt[f38,μm[f38]]+#,μm[f38],0,True,(10^(-4)),False,4]&/@Range[0,25/1000,5/1000];

testMaxMarge=solLineMin[f38,Eopt[f38,μm[f38]]+#,μm[f38],0,False,(10^(-4)),False,4]&/@Range[0,25/1000,5/1000];

testMaxμMarge=solLineMinμ[f38,Eopt[f38,μm[f38]]+#,μm[f38]-1/10000,0,True,(10^(-2)),False,3]&/@Range[0,25/1000,5/1000];

testMinμMarge=solLineMinμ[f38,Eopt[f38,μm[f38]]+#,μm[f38]-1/10000,0,False,(10^(-2)),False,3]&/@Range[0,25/1000,5/1000];

testDomMargMarg=solLineDominantFixedμ[f38,Eopt[f38,μm[f38]]+#,μm[f38],μm[f38],(10^(-4))]&/@Range[5/1000,25/1000,5/1000];

testDomMarg=solLineDominant[f38,Eopt[f38,μm[f38]]+#,μm[f38],(10^(-3))]&/@Range[0,25/1000,5/1000];

ListLogLogPlot[{1-q,E1-E0}/.testDomMargMarg,Joined->True]

test5=addEigenvalueInformation/@testMinMarge2;

test4=addEigenvalueInformation/@testDomMargMarg;

ListLogLogPlot[{1-q,E1-E0}/.(Select[#,μ1-(y+(1/y)f''[1])<0&&y>μm[f]/2/.#&]&/@test4),Joined->True]

p4=Show[ListLogLogPlot[{1-q,E1-E0}/.MapThread[Join[Most@#1,Reverse[Most@#2]]&,{testMinMarge2,testMaxMarge2}],Joined->True,Frame->True,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[E, 1])"}],RotateLabel->False]

eigMinMarge=addEigenvalueInformation/@testMinMarge;

p4=Show[ListLogLogPlot[{1-q,E1-E0}/.MapThread[Join[Most@#1,Reverse[Most@#2]]&,{testMinMarge,testMaxMarge}],Frame->True,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[E, 1])"}],RotateLabel->False]

Sin[ArcCos[q]/2]//FullSimplify[#,Assumptions->{q>0,q<1}]&

2n-2 q n

rtest1=Log[((π^(n/2))/(Γ[n/2+1]))(r^n)]/.r->(Sqrt[2n])(Sqrt[1-q])/2/.Γ[x_]:>(Sqrt[2π])(x^(x-1/2))Exp[-x]//FullSimplify[#,Assumptions->{n>2,q<1,q>0}]&

rtest2=Limit[rtest1/n,n->∞]

ListPlot[{E0,q}/.Select[#,y<μm[f]/2/.#&][[-1]]&/@test5]

ListPlot[rtest2/.Select[#,y<μm[f]/2/.#&][[-1]]&/@test5]

EE/.groundState[f38,EE,μμ]

testMinMarge2[[1,1]]//N

ListPlot[{((E0-(EE/.groundState[f38,EE,μμ]))/(Abs[EE/.groundState[f38,EE,μμ]])),rtest2+(𝒟[f][μ0]+𝒮1[f][E0,μ0])}/.Select[#,y<μm[f]/2/.#&][[-1]]&/@test5]

p4=Show[ListLogLogPlot[{1-q,E1-E0}/.testDomMarg,Joined->True,Frame->True,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[E, 1])"}],RotateLabel->False]

p4=Show[ListLogLogPlot[{1-q,Abs[E1-E0]}/.Rest@testMaxMarge,Joined->True,Frame->True,Prolog->{{Dashed,Line[{{0.8,Eopt[f34,μm[f34]]+#},{1,Eopt[f34,μm[f34]]+#}}]&/@Delete[Rest@Range[-2/100,1/100,5/1000],{4}]},Dashed,Line[{{0.4,Eopt[f34,μm[f34]]},{1,Eopt[f34,μm[f34]]}}],{Directive[FontFamily->Times],Text["(Subscript[E, th])",{0.35,Eopt[f34,μm[f34]]}]},{Directive[FontFamily->Times],Text["(Subscript[E, 0])",{0.9,Eopt[f34,μm[f34]]-18/1000}]}},ImageSize->250,AspectRatio->1,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[E, 1])"}],RotateLabel->False]

p5=Show[ListPlot[{q,μ1}/.testMinMarge,Joined->True,Frame->True,Prolog->{Dashed,Line[{{0.4,μm[f38]},{1,μm[f38]}}],{Directive[FontFamily->Times],Text["(Subscript[μ, m])",{0.35,μm[f34]}]},{Directive[FontFamily->Times],Text["(Subscript[E, 0])",{0.9,Eopt[f34,μm[f34]]-18/1000}]}},ImageSize->700,AspectRatio->1,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[μ, 1])"}],RotateLabel->False]

test1=addEigenvalueInformation/@testMaxμMarge;

f38''[1]

μm[f38]//N

ListLogPlot[{q,y+(1/y)f38''[1]-μm[f38]}/.(Select[#,y>μm[f38]/2/.#&]&/@test1),Joined->True]

((E1-E0)/E1)/.gs38/.test1[[All,1]]

Show[
ListLogLogPlot[{1-q,Abs[6-μ1]}/.(testMinMarge),Joined->True,Frame->True,PlotRange->{{0,0.1},{-(10^(-2)),(10^(-2))}},ImageSize->700,AspectRatio->1,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[μ, 1])"}],
RotateLabel->False,Prolog->Line[{{0,μm[f38]},{1,μm[f38]}}]
]

Show[ListLogPlot[{1-q,E1-E0}/.MapThread[Join[Most@#1,Reverse[Most@#2]]&,{testMaxμMarge,testMinμMarge}],Joined->True,Frame->True,PlotRange->{{0,0.1},{-(10^(-2)),(10^(-2))}},Prolog->{Dashed,Line[{{0.4,μm[f38]},{1,μm[f38]}}],{Directive[FontFamily->Times],Text["(Subscript[μ, m])",{0.35,μm[f38]}]}},ImageSize->700,AspectRatio->1,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[μ, 1])"}],RotateLabel->False,Prolog->Line[{{0,μm[f38]},{1,μm[f38]}}]]

f38

Show[ListPlot[{q,μ1}/.testMinμMarge,Joined->True,Frame->True,PlotStyle->(Blend[{ColorData[97][2],ColorData[97][3]},#]&/@Range[0,1,1/12]),Prolog->{Dashed,Line[{{0.4,μm[f34]},{1,μm[f34]}}],{Directive[FontFamily->Times],Text["(Subscript[μ, m])",{0.35,μm[f34]}]},{Directive[FontFamily->Times],Text["(Subscript[E, 0])",{0.9,Eopt[f34,μm[f34]]-18/1000}]}},ImageSize->700,AspectRatio->1,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[μ, 1])"}],RotateLabel->False]

ColorData[97]/@Range[10]

p5=Show[ListLogPlot[{q,μ1-μm[f38]}/.testMinMarge,Joined->True,Frame->True,ImageSize->700,AspectRatio->1,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[μ, 1])"}],RotateLabel->False]

Show[ListLogPlot[{q,μ1-μm[f38]}/.testMaxμMarge,Joined->True,Frame->True,PlotStyle->(Blend[{ColorData[97][2],ColorData[97][3]},#]&/@Range[0,1,1/10]),ImageSize->700,AspectRatio->1,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[μ, 1])"}],RotateLabel->False]

Show[ListLogLogPlot[{1-q,Abs[μm[f38]-μ1]}/.testMinμMarge,Joined->True,Frame->True,ImageSize->700,AspectRatio->1,LabelStyle->Directive[Black,FontFamily->Times],FrameLabel->{q,"(Subscript[μ, 1])"}],RotateLabel->False]

ListLogLogPlot[{1-q,Abs[μ1-μ0]}/.testMinMarge,Joined->True]

Evaluate[D[((f38''[q])^(-1/2)),{q,2}]]

FindMinimum[D[((f38''[q])^(-1/2)),{q,2}],{q,0.6},WorkingPrecision->40]

gs38=groundState[f38,E1,μ1]

0.01 E1/.gs38

(E1-Eopt[f38,μm[f38]])/E1/.gs38

Plot[((f38''[q])^(-1/2)),{q,0,1}]

qMCT38=FindRoot[1-q==((f38'[q])/(q f38''[q])),{q,0.5},WorkingPrecision->30]

TMCT38=(Sqrt[((1-q)/q)f38'[q]])/.qMCT38

qMCT34=FindRoot[1-q==((f34'[q])/(q f34''[q])),{q,0.5},WorkingPrecision->30]

TMCT34=(Sqrt[((1-q)/q)f34'[q]])/.qMCT34

LogPlot[Evaluate[D[((f38''[q])^(-1/2)),{q,2}]],{q,0,1},PlotRange->All]

dat=Import["/tmp/Mixed_spin_p3q8_lambda0.76_Tf0_Ti0.542345_L5000_dt0.05.dat"];

fitμ=NonlinearModelFit[{6-#6,#4}&@@@dat,a+b ((Δμ)^c),{a,b,c},Δμ]

fitE=NonlinearModelFit[Rest[{#1,#4}&@@@dat],a+b (t^(-c)),{a,b,c},t]

Show[
ListPlot[{6-#6,#4}&@@@dat,PlotRange->All],
Plot[fitμ[x],{x,0,2}]
]

Show[
ListLogLinearPlot[{#1,#4}&@@@dat,PlotRange->All],
LogLinearPlot[fitE[x],{x,0,10000}]
]

underRange=Range[-15/1000,-25/10000,25/10000];

Eopt[f34,μm[f34]]+Range[-15/1000,-25/10000,25/10000]//N

testMinMargeBelow=solLineMin[f34,Eopt[f34,μm[f34]]+#,μm[f34],0,True,(10^(-3))]&/@underRange;

testMaxMargeBelow=solLineMin[f34,Eopt[f34,μm[f34]]+#,μm[f34],0,False,(10^(-3))]&/@underRange;

testMinMarge2Below=solLineMarginalMin[f34,Eopt[f34,μm[f34]]+#,μm[f34],0,True,(10^(-3))]&/@underRange;

testMaxMarge2Below=solLineMarginalMin[f34,Eopt[f34,μm[f34]]+#,μm[f34],0,False,(10^(-3))]&/@underRange;

testMaxμMargeBelow=solLineMinμ[f34,Eopt[f34,μm[f34]]+#,μm[f34],0,True,(10^(-4))]&/@underRange;

testMinμMargeBelow=solLineMinμ[f34,Eopt[f34,μm[f34]]+#,μm[f34],0,False,(10^(-3))]&/@underRange;

testDomMargMargBelow=solLineDominantFixedμ[f34,Eopt[f34,μm[f34]]+#,μm[f34],μm[f34],(10^(-3))]&/@underRange;

testDomMargBelow=solLineDominant[f34,Eopt[f34,μm[f34]]+#,μm[f34],(10^(-3))]&/@underRange;

ListPlot[{1-q,E1-Eopt[f,μm[f]]}/.testDomMargMargBelow,Joined->True]

test2=addEigenvalueInformation/@testDomMargMargBelow;

test3=addEigenvalueInformation/@testMaxMarge2;

test2[[4,-5]]

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

(Eopt[f,μm[f]]/.f->f34)+underRange//N

test2[[1,1]]//N

Esf=((-169981)/100000);

testDomSF=solLineDominantFixedμ[f34,Esf,μm[f34],μm[f34],(10^(-3))];

testEigSF=addEigenvalueInformation@testDomSF;

testEigMinMarge2=addEigenvalueInformation/@testMinMarge2;

ListPlot[{1-q,E1-Eopt[f,μm[f]]}/.(Select[#,μ1-(y+(1/y)f''[1])<0&&y>μm[f]/2/.#&]&/@MapThread[Join,{test2,Reverse/@test3}]),Joined->True,PlotStyle->Dashed]

Show[ListLogPlot[{1-q,E1-E0}/.(Select[#,μ1-(y+(1/y)f''[1])<0&&y>μm[f]/2/.#&]&/@MapThread[Join,{test3,Reverse/@testEigMinMarge2}]),Joined->True,PlotStyle->Dashed],
ListLogPlot[{1-q,E1-E0}/.(Select[#,y<μm[f]/2/.#&]&/@MapThread[Join,{test3,ConstantArray[{{y->0,E1->I}},Length[test3]],Reverse/@testEigMinMarge2}])]
]

Show[ListLogPlot[{1-q,E1-E0}/.(Select[#,y<μm[f]/2/.#&]&/@test2),Joined->True],
ListLogPlot[{1-q,E1-E0}/.(Select[#,μ1-(y+(1/y)f''[1])<0&&y>μm[f]/2/.#&]&/@test2),Joined->True,PlotStyle->Dashed]]

test=(Select[#,μ1-(y+(1/y)f''[1])<0&&y>μm[f]/2/.#&]&@testEigSF)[[1]]

Eopt[f34,μm[f34]]//N

dettest=Det[Simplify[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]/.{y->μm[f]/2,μ0->μm[f],μ1->μm[f],E1->Eopt[f,μm[f]],β1->0}/.f->f34,Assumptions->{f''[1]>0,f'[1]>0,f''[1]>f'[1]}]]

SFdat=FindRoot[Evaluate@Append[D[𝒮12[f34][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{{β1,r11d,r110,r01,r10,q0}}]/.{β1->0,μ0->μm[f],μ1->μm[f],E1->Eopt[f,μm[f]]}/.f->f34,dettest],test[[{2,5,8,9,10,11,12}]]/.Rule->List,WorkingPrecision->30]

dettest/.SFdat

D[𝒮12[f34][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0],{{β1,r11d,r110,r01,r10,q0}}]/.{β1->0,μ0->μm[f],μ1->μm[f],E1->Eopt[f,μm[f]]}/.f->f34/.SFdat

ESF[f34]//N

E0-ESF[f34]/.SFdat

ListPlot[{1-q,E1-Eopt[f,μm[f]]}/.(Select[#,μ1-(y+(1/y)f''[1])<0&&y>μm[f]/2/.#&]&@testEigSF),Joined->True,PlotStyle->Dashed]

Show[ListPlot[{1-q,E1-Eopt[f,μm[f]]}/.(Select[#,μ1-(y+(1/y)f''[1])<0&&y>μm[f]/2/.#&]&/@MapThread[Join,{test3,Reverse/@testEigMinMarge2}]),Joined->True,PlotStyle->Dashed],
ListLogPlot[{1-q,E1-Eopt[f,μm[f]]}/.(Select[#,y<μm[f]/2/.#&]&/@MapThread[Join,{test3,ConstantArray[{{y->0,E1->I}},Length[test3]],Reverse/@testEigMinMarge2}])]
]

ListLogPlot[{1-q,E1-E0}/.MapThread[Join,{testMinMargeBelow,Reverse/@testMaxMargeBelow}],Joined->True,PlotRange->{-0.01,0.01}]

ListPlot[{1-q,E1-Eopt[f,μm[f]]}/.(Select[#,y>μm[f34]/2&&𝒟[f][μ1]+𝒮12[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0]>0/.#&]&/@test2),Joined->True]

Show[
ListLogLogPlot[{1-q,E1-E0}/.(Select[#,y<μm[f34]/2/.#&]&/@test2),Joined->True],
ListLogLogPlot[{1-q,E1-E0}/.(Select[#,y>μm[f34]/2/.#&]&/@test2),Joined->True,PlotStyle->Dashed]
]

test6=addEigenvalueInformation/@testDomMargBelow;

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

ListLogLogPlot[{1-q,Abs[E1-Eopt[f,μm[f]]]}/.(Select[#,𝒟[f][μ1]+𝒮12[f][E0,μ0,E1,μ1,q][β1,r11d,r110,r01,r10,q0]>0&&μ1>μ0&&μ1-(y+(1/y)f''[1])<0&&y>μm[f]/2/.#&]&)/@test6,Joined->True]

1

Show[
ListLogLogPlot[{1-q,E1-E0}/.(Select[#,y+(1/y)f38''[1]-μ1<0/.#&]&/@test3),Joined->True,PlotStyle->Dashed],
ListLogLogPlot[{1-q,E1-E0}/.(Select[#,y<μm[f38]/2/.#&]&/@test3),Joined->True]
]